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ABSTRACT 

This  report  serves  to  document  progress  made  to  date  on  National  Science  Foundation 
Project  #97-13481,  Earthquake  Hazard  Mitigation  Program.  The  focus  of  this  phase  of  the 
project  is  the  development  of  an  improved  solution  algorithm  for  fast  transient  analysis  of  large, 
locally  nonlinear  structures  using  time  domain  structural  synthesis.  Time  domain  structural 
synthesis  is  a  general  and  exact  formulation  for  transient  problems  in  structural  modification, 
substructure  coupling,  and  base  excitation.  The  formulation  is  characterized  by  the  governing 
equation  of  the  synthesis,  which  is  a  nonlinear  Volterra  integral  equation.  The  governing 
equation  makes  use  of  impulse  response  functions  calculated  for  those  coordinates  of  the 
(sub)structures  subjected  to  forces  of  synthesis  (e.g.  modification  forces,  coupling  forces).  This 
physical  coordinate  formulation  provides  for  a  largely  unrestricted  and  exact  model  reduction,  in 
that  only  those  coordinates  of  interest  need  be  retained  in  the  synthesis. 

The  report  documents  the  development  of  several  algorithms  intended  to  improve  upon 
the  original  algorithm  developed  by  the  first  author.  The  last  algorithm  developed  is  presented 
first  in  this  report,  as  this  algorithm  meets  the  stated  goals  of  the  project.  This  algorithm  is  a 
newly  developed  recursive,  block-by-block  convolution  solution  to  the  governing  nonlinear 
integral  equation.  As  is  demonstrated  with  a  simple  but  realistically  large  nonlinear  base 
excitation  problem  (51,500  DOF  finite  element  model),  the  new  algorithm  provides  a  78% 
reduction  in  time  required  for  the  nonlinear  transient  base  excitation  solution,  as  compared  with 
traditional  direct  integration  calculated  using  a  widely-used  commercial  finite  element  program. 
This  very  large  savings  in  computer  time  is  obtained  for  a  single  analysis,  i.e.  assuming  no  prior 
calculations  have  been  made  for  the  impulse  response  functions  of  the  (sub)structures.  The  new 
algorithm  provides  an  even  greater  reduction  in  computer  time  for  all  subsequent  analyses.  As 
shown  in  the  example  problem,  once  all  required  impulse  response  functions  have  been 
calculated,  the  nonlinear  base  isolation  solutions  calculated  using  the  new  recursive,  block-by- 
'  block  convolution  algorithm  take  approximately  7  seconds,  as  compared  with  the  direct 
integration  solution  which  takes  approximately  30  minutes.  This  rapid  reanalysis  capability  will 
facilitate  the  development  of  numerical  optimization  for  the  design  of  nonlinear  isolation.  The 
theory  of  transient  synthesis  is  documented,  along  with  a  new  proof  of  the  exponential 
convergence  properties  of  an  iterative  solution  to  the  governing  nonlinear  integral  equation. 
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1. 


INTRODUCTION 


The  transient  analysis  of  large  and  complex  structural  systems  is  a  computationally 
demanding  task  exacerbated  by  the  presence  of  structural  and  mechanical  nonlinearities.  The 
computational  demand  of  these  problems  prohibits  the  repeated  analyses  required  in  a  design 
effort,  such  as  in  structural  optimization  where  various  responses  are  required  for  the  calculation 
of  the  objective  function,  constraints,  sensitivities,  and  for  the  generation  of  approximations  to  be 
used  within  the  optimizer. 

A  class  of  nonlinear  structural  dynamics  problems  with  numerous  applications  is 
characterized  by  the  presence  of  localized,  nonlinearities.  For  the  purposes  of  this  work,  this  class 
of  problems  is  defined  as  follows: 

Definition  of  a  Locally  Nonlinear  Model:  A  model  where  the  nonlinear  load  paths 

do  not  contain  any  internal  degrees-of-freedom  (DOF),  i.e.  each  nonlinear  load  path 
(nonlinear  element)  is  associated  solely  with  DOF  shared  by  linear  load  paths  (elements). 

This  class  of  problems  can  be  further  informally  restricted  by  recognizing  that  the  formulations 
to  be  developed  in  what  follows  provide  a  greater  reduction  in  computing  time  (as  compared 
with  direct  integration)  for  models  where  there  are  relatively  few  nonlinear  load  paths,  or  in  other 
words,  where  the  number  of  DOF  associated  with  nonlinear  load  paths  is  small  relative  to  the 
total  number  of  DOF  in  the  model.  The  problem  of  nonlinear  earthquake  isolation  of  a  linear 
structure  falls  into  this  category,  wherein  the  isolator  provides  a  nonlinear  load  path  between  the 
building  model  DOF  and  “ground.” 

A  locally  nonlinear  model  allows  for  algorithmic  approaches  which  exploit  this 
characteristic  to  achieve  reduced  computer  time  requirements.  The  most  common  strategy  for 
this  exploitation  is  where  the  localized  nature  of  the  nonlinearities  facilitates  their  treatment  as 
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terms  on  the  right  hand  side,  i.e.  applied  loads  of  differential  equations,  albeit  loads  which  are 
functions  of  system  responses,  and  possibly  with  time,  see  for  example  [1]. . 

The  approach  taken  in  this  work  is  to  treat  the  problem  as  a  physical  coordinate  (non- 
modal)  structural  modification  problem,  wherein  the  nonlinear  elements  are  “installed”  into  the 
linear  model  as  structural  modifications.  The  structural  modification  formulation  belongs  to  the 
broader  category  of  physical  coordinate  structural  synthesis  [2-6],  which  includes  substructure 
coupling  and  constraint  imposition  as  well.  Such  an  approach  not  only  provides  a  substantial 
reduction  in  solution  times,  but  provides  for  a  generality  in  the  definition  of  the  problem  and  a 
flexibility  in  its  application  which  is  unique.  We  will  describe  structural  synthesis  here. 

In  a  sense,  structural  synthesis  treats  the  nonlinear  element  responses  as  applied  loads  as 
well.  However,  what  distinguishes  structural  synthesis  from  other  numerical  approaches  are  the 
following  characteristics: 

•  The  governing  equations  for  structural  synthesis  are  exact. 

•  An  implicit  exact  model  reduction  is  available,  in  that,  as  a  minimum  only  those  DOF 

directly  associated  with  nonlinear  elements  and  applied  loading  need  be  retained.  Any 
additional  physical  DOF  of  interest  to  the  analyst  can  be  retained  as  well.  In  other  words, 
the  transient  synthesis  solution  time  is  independent  of  model  size. 

•  Very  general  nonlinearities  can  be  treated. 

•  The  linear  portion(s)  of  the  model  is  solved  once. 

•  Very  fast  solution  times  are  obtained,  an  intrinsic  property  of  the  formulation. 

1.1  Summary  of  Results  Reported  for  Year  One 

The  goal  of  the  Year  One  phase  of  this  project  is  stated  here: 

The  development  of  a  highly  efficient  and  general  formulation  for  nonlinear  transient 
structural  synthesis  which  will  provide  fast  linear  or  nonlinear  transient  re-analysis  with 
arbitrary  loading  for  large  linear  structural  FE  models  which  can  have  localized,  but 
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generally  nonlinear  components  of  arbitrary  magnitude.  The  nonlinear  formulation  will  be  an 
extension  of  prior  research  of  the  PI  [3]. 

This  goal  is  met  by  the  algorithm  developed  in  Section  3  of  this  report. 

1.2  Notes  on  Algorithms  Developed  -  Organization  of  Report 

The  work  reported  herein  results  in  the  development  of  a  fast  numerical  solution 
algorithm  for  the  governing  equation  of  locally  nonlinear  transient  structural  synthesis.  This  is 
the  primary  goal  of  the  Year  One  effort.  In  pursuit  of  this  goal,  three  general  recursive  algorithms 
were  developed,  and  these  algorithms  are  numbered  in  the  order  in  which  they  were  developed. 
However,  the  third  algorithm  developed  resulted  in  the  recursive  block-by-block  convolution 
algorithm  (Recursive  Algorithm  3)  which  meets  the  stated  goals  (Section  1.1)  for  this  phase  of 
the  project.  This  algorithm  will  be  the  focus  of  continued  work,  and  this  algorithm  will  be 
presented  first  (Section  3). 

We  will  detail  all  of  the  algorithms  developed,  including  those  either  deemed  ineffective 
for  the  task  at  hand  (nonlinear  isolation  simulation)  or  surpassed  by  algorithms  subsequently 
developed  (and  described  herein)  for  the  purpose  of  documenting  “lessons  learned,”  and  to 
document  novel  results  obtained,  of  potential  value  in  other  areas.  We  summarize  briefly  the 
algorithms  here. 

The  governing  equation  for  transient  structural  synthesis  is  a  nonlinear  Volterra  integral 
equation,  involving  a  convolution-type  kernel  [3],  The  convolution-type  kernel  suggests  the 
recursive  transition-matrix  approach  to  the  solution  of  first-order  ordinary  differential  equations 
as  a  potential  improvement  Over  the  (non-recursive)  iteration  solution  presented  by  Gordis  and 
Radwick  [3],  which  itself  demonstrated  an  order-of-magnitude  reduction  in  computing  time 
required,  relative  to  direct  nonlinear  transient  integration.  A  recursive  transition-matrix  approach 
is  taken  in  the  work  of  Inaudi  and  De  La  Llera  [7],  and  this  work  provided  additional  motivation 
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in  this  direction.  However,  in  [7],  the  recursion  was  based  on  the  transition  matrix  for  the  system 
model,  and  hence  requires  the  calculation  of  a  large  matrix  exponential,  with  no  provision  for 
model  reduction.  We  must  therefore  consider  such  an  approach  to  be  of  limited  value  for  large 
structural  models.  Our  goal  was  to  explore  the  development  of  a  recursive  transition-matrix 
algorithm  appropriate  for  large  models,  i.e.  having  an  implicit  model  reduction.  The 
developments  and  results  of  this  phase  of  the  work  are  detailed  in  Section  4,  and  include 
Recursive  Algorithms  1  and  2,  which  are  summarized  here. 

The  first  two  of  the  above-mentioned  algorithms  are  based  on  this  transition-matrix 
approach.  These  algorithms  are  as  follows: 

1.2. 1  Recursive  Algorithm  1  (RA 1 ) 

RA1  is  included  as  a  necessary  precursor  to  the  algorithm  RA2  which  follows,  rather  than 
due  to  any  particular  novelty  or  computational  value.  The  algorithm  RA1  is  based  on  transition 
matrices  constructed  from  the  modal  transient  response  in  first-order  (state-space)  form.  As  will 
be  shown,  this  algorithm  is  inherently  unstable,  and  can  only  be  made  stable  by  the  inclusion  of 
unrealistic  levels  of  system  damping.  It  is  included,  however,  as  a  precursor  for  the  following 
algorithm,  and  as  a  baseline  for  comparison  purposes. 

1.2.2  Recursive  Algorithm  2  (RA2) 

This  algorithm  develops  a  new  recursive  transition  matrix  solution  for  transient  structural 
dynamics.  The  algorithm  defines  a  complex  modal  state  vector,  and  an  associated  transition 
matrix  which  is  based  on  second-order  modal  differential  equations.  The  motivation  for  the 
development  of  this  algorithm  was  to  avoid  the  formulation  of  the  problem  in  terms  of  first-order 
differential  equations,  i.e.  such  as  in  RA1,  in  order  to  maintain  consistency  of  the  formulation 
with  common  tools  for  structural  dynamics  analysis,  e.g.  commercial  finite  element  programs. 
Furthermore,  as  will  be  shown,  RA2  achieves  a  significant  computational  advantage  relative  to 
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RA1,  as  RA2  is  more  “diagonalized.”  Unfortunately,  RA2  is  also  inherently  unstable,  again,  only 
stabilized  by  the  inclusion  of  substantial  system  damping. 

1.2.3  Recursive  Predictor-Corrector  Algorithms  1  and  2  (RA1PC  and  RA2PC) 

The  inherent  instability  of  algorithms  RA1  and  RA2  motivated  their  reformulation  as  predictor- 
corrector  algorithms,  similar  to  the  approach  in  [7].  We  will  refer  to  these  implicit  versions  as 
RA1PC  and  RA2PC. 

Algorithms  RA1  and  RA2  are  reformulated  as  predictor-corrector  algorithms,  with  a 
single  corrector  step.  The  stability  calculations  for  RA1PC  and  RA2PC  indicate  that  it  is  possible 
to  stabilize  these  algorithms  by  going  to  implicit  versions.  However,  the  results  indicate  that 
further  development  of  these  algorithms  would  be  necessary  in  order  to  improve  on  the 
performance  already  obtained  by  Gordis  and  Radwick  [3].  Hence,  Recursive  Algorithm  3  is 
developed. 

1.2.4  Recursive  Algorithm  3:  Block-by-Block  Convolution 

This  algorithm  is  markedly  different  from  RA1  and  RA2  in  that  no  transition  matrix  is 
employed.  This  algorithm  preserves  the  physical  coordinate  formulation  originally  developed  by 
Gordis  [2]  and  Gordis  and  Radwick  [3],  and  hence  preserves  the  implicit  and  unrestricted  exact 
model  reduction,  concomitant  with  the  formulation.  The  algorithm  will  be  shown  to  be 
exponentially  convergent,  for  a  general  class  of  nonlinearities.  The  current  version  of  the 
software  makes  use  of  an  FIR  filter  to  calculate  the  convolution  sum,  as  implemented  in  the 
MATLAB®  built-in  convolution  function  (“CONV”).  It  should  be  noted  that  the  FIR 

implementation  in  MATLAB®  is  not  optimally  efficient  as  one  can  develop  a  filter  which 
computers  only  those  elements  required.  However,  this  characteristic  of  MATLAB ’s  convolution 
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algorithm  will  yield  conservative  computer  time  comparisons.  Recursive  Algorithm  3  (RA3) 
the  first  algorithm  discussed,  in  Section  3  of  this  report. 


2. 


INTEGRAL  EQUATION  FORMULATION  FOR 
TRANSIENT  STRUCTURAL  SYNTHESIS 


As  the  Year  One  task  focuses  on  the  development  of  an  improved  solution  algorithm  for 
the  governing  equation  for  locally  nonlinear  transient  structural  synthesis,  we  provide  the 
background  in  the  relevant  theory.  The  reader  is  referred  to  Gordis  [2]  and  Gordis  and  Radwick 
[3]  for  the  complete  development.  The  following  section  is  abstracted  from  [3]. 

The  theory  defines  a  transient  analysis  that  is  independent  of  model  size,  in  that  the 
theory  is  cast  in  physical  coordinates  (i.e.  non-modal),  and  the  transient  analysis  is  done  using 
only  those  structural  DOF  of  interest.  These  DOF  must  include,  as  a  minimum,  those  associated 
with  the  nonlinear  elements,  which  are  treated  independently  of  the  (linear)  model.  Additionally, 
other  DOF  for  which  synthesized  response  information  is  desired  can  be  included  as  needed. 
Therefore,  it  is  possible  to  synthesize  the  transient  response  for  an  arbitrarily  large  model  using  a 
minimal  number  of  DOF,  the  minimum  number  defined  only  by  the  number  of  nonlinear 
elements  in  the  model.  This  unique  feature  of  the  theory  has  been  demonstrated  for  the  linear 
transient  formulation  [2],  the  nonlinear  transient  formulation  [3],  and  in  the  frequency  domain  in 
[4-6].  Functioning  as  a  re-analysis  procedure,  the  formulation  directly  calculates  the  new 
transient  response  for  a  system  resulting  from  structural  changes  and/or  coupling  with  other 
structures,  without  a  reassembly  or  full  reanalysis.  These  features  will  be  more  fully  described  in 
what  follows. 

Structural  synthesis  also  provides  for  substructure  coupling.  This  feature  allows  nonlinear 
elements  to  be  isolated  by  division  of  the  system  into  substructures,  which  provides  additional 
computational  advantages  in  that  the  synthesis  can  exploit  inherent  physical  boundaries  in  a 
problem.  One  immediate  benefit  resulting  is  that  linear  substructures  are  solved  once,  which 
constitutes  a  significant  savings  in  that  the  eigensolution  is  typically  the  most  expensive  part  of  a 
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dynamic  analysis  [8].  The  synthesis  is  used  to  connect  the  substructures  through  the  nonlinear 
elements,  and  to  calculate  the  combined  system  nonlinear  transient  response. 

Each  substructure  is  described  by  impulse  response  functions  calculated  at  the  DOF 
where  nonlinear  elements  are  to  be  installed,  where  loads  are  applied,  and  at  other  DOF  for 
which  synthesized  nonlinear  transient  response  is  required.  For  each  linear  substructure,  the  IRF 
are  most  efficiently  calculated  using  modal  superposition.  However,  the  use  of  modal 
superposition  for  IRF  calculation  does  not  render  structural  synthesis  a  “modal  method,”  for  the 
following  reason.  The  IRF  are  calculated  using  a  sufficient  number  of  modes  to  ensure 
convergence.  Once  these  converged  IRF  are  calculated,  they  are  indistinguishable  (to  a  given 
level  of  precision)  from  the  “exact”  IRF,  which  are  indeed  physical  quantities.  We  can  contrast 
this  approach  with  modal  methods  of  structural  synthesis,  where  the  convergence  relative  to  the 
number  of  retained  modes  is  determined  from  the  solution  of  the  synthesized  system.  Typically, 
a  modal  synthesis  procedure  generates,  from  substructure  solutions,  a  coupled  system  model  for 
which  an  additional  solution  is  required.  Furthermore,  modal  synthesis  methods  are  inherently 
approximate.  The  governing  equations  of  structural  synthesis  to  be  developed  here  are  exact.  The 
only  errors  incurred  are  those  avoidable  errors  resulting  from  retaining  an  insufficient  number  of 
modes,  and  those  errors  common  to  all  numerical  integration  schemes,  which  are  discretization 
errors,  and  truncation  errors  where  numerical  derivatives  are  used. 

2.1  Notation  Used 

In  what  follows,  all  vector  quantities  are  denoted  by  boldface  type,  e.g.  x.  A  scalar 
element  of  a  vector  quantity  will  be  denoted  by  a  subscripted  plainface  type,  e.g.  xi5  which 
denotes  the  i*  element  of  the  vector  x.  We  will  define  various  coordinate  sets,  which  often  are 
defined  as  subsets  of  a  vector  of  coordinates.  For  example,  if  the  coordinate  vector  x  is 
comprised  of  two  subsets  of  coordinates,  say  the  a  set  and  the  (3  set,  we  will  denote  these 
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coordinate  subsets  as  xa  and  xp.  The  coordinate  set  contained  in  the  vector  x  is  the  union  of 
subsets  a  and  (3.  The  vector  x  can  be  written  as  x  =  [^xj  xp  J  .  A  particular  vector  is  denoted 
by  a  subscripted  boldface  symbol,  such  as  the  i*  eigenvector,  denoted  <jv 

We  will  also  need  to  distinguish  quantities  associated  with  various  mathematical  domains 
such  as  physical,  modal,  state-space,  and  the  vector  quantities  comprised  thereof.  Quantities 
comprised  of,  or  associated  with  physical  (nodal)  coordinates  will  have  no  embellishment. 
Quantities  comprised  of,  or  associated  with  modal  coordinates  will  be  denoted  with  the  tilde 
embellishment  (•).  Quantities  comprised  of,  or  associated  with  state-space  coordinates  will  be 
denoted  with  the  caret  embellishment  (•).  Quantities  comprised  of,  or  associated  with  the 
complex  modal  coordinates  (to  be  defined)  will  be  denoted  with  the  bar  embellishment  (•) .  In 
general,  these  will  also  be  evident  from  context. 

The  various  algorithms  presented  involve  discrete  time  histories,  e.g.  x(kAt)  or  x(kT) 

where  the  sample  length  is  At  =T  in  seconds.  We  will  abbreviate  this  notation  with  the  following 
superscript  notation,  x(kAt)  =  x(kT)  =  x00. 

2.2  Theory 

In  the  context  of  the  physical  coordinate  synthesis  formulation  to  be  developed,  a 
structural  system  is  defined  to  consist  of  one  or  more  uncoupled  substructures.  A  single 
governing  equation  for  nonlinear  transient  synthesis  will  be  derived  and  this  equation  will 
address  each  of  the  following  three  general  analysis  categories: 

(1)  Structural  modification  -  the  addition  and/or  removal  of  linear  and/or  nonlinear  structural 
elements 

(2)  Prescribed  base  motion  -  application  of  base  motion  to  structure  through  linear  and/or 
nonlinear  elements 
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(3)  Substructure  coupling  -  the  joining  of  substructures  (a  linear  analysis) 


Each  of  the  above  analysis  categories  defines  a  set  of  DOF.  Referring  to  Figure  1,  a  structural 
system  is  shown,  comprised  of  two  substructures,  each  of  an  arbitrary  number  of  DOF. 


The  subset  of  DOF  associated  with  structural  modifications  is  denoted  the  “m-set”  and 
may  include  DOF  from  all  substructures.  The  subset  of  DOF  associated  with  prescribed  base 
motion  excitation  is  denoted  the  “b-set”  and  again,  may  include  DOF  from  all  substructures.  The 
subset  of  DOF  associated  with  substructure  coupling  is  denoted  the  “c-set”  and  may  include 
DOF  from  all  substructures.  Given  that  all  nonlinear  elements  are  installed  in  the  synthesis,  the 
substructures  are  linear,  and  hence  coupling  is  a  linear  synthesis  [2].  The  subset  of  DOF  denoted 
the  “i-set”  refers  to  those  system  DOF  about  which  synthesized  response  information  is  required, 
but  are  not  directly  involved  with  the  synthesis,  either  modification,  base  excitation,  or  coupling. 
The  DOF  sets  are: 

m-set:  Modification 
b-set:  Base  excitation 
c-set:  Coupling 
i-set:  Additional  DOF 
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While  much  of  the  attention  in  this  work  will  focus  on  single  substructure  systems,  we 
will  often  refer  to  “substructures”  as  opposed  to  “structures”  in  keeping  with  the  formulation  of 
the  problem  as  one  of  structural  synthesis. 

2.3  Basic  Definitions  for  a  Substructure  Model 

A  finite  element  (FE)  discretization  produces  the  following  linear  (sub)structure  model: 

Mx  +  Cx  +  Kx  =  fe(t)  (2.1) 

where  the  mass  matrix  M  is  (N  x  N)  symmetric  positive  definite,  the  damping  matrix  C  is  (N  x 
N)  symmetric  positive-semidefinite,  and  the  stiffness  matrix  K  is  (N  x  N)  symmetric  positive- 
semidefinite.  The  inertial  nodal  displacement  vector  x(t)  is  (N  x  1),  as  are  the  nodal  velocity  and 
acceleration,  denoted  x(t)  and  x(t)  respectively.  The  loading  vector  fe(t)  contains  time- 
dependent  loads  acting  at  the  nodal  degrees-of-freedom  (DOF).  The  FE  model  possesses  N  DOF. 
In  the  case  of  base  excitation  (e.g.  ground  motions),  the  loading  vector  can  be  written  as, 

4(t)  =  Kgy(t)  +  Cgy(t)  (2.2) 

where  Kg  and  Cg  are  the  stiffness  and  damping  matrices  associated  with  (multiple)  base 
displacement  and  velocity  excitation.  Note  that  applied  forces  and  moments,  and  base  excitation 
may  co-exist. 

The  mass  and  stiffness  matrix  possess  the  eigenpairs,  (Oj,^,  i=l,2,..,N.  The  eigenvalue 
a)  has  the  units  sec'1.  The  structure  may  possess  up  to  6  rigid-body  modes,  distinguished  by  co—0, 
i  =  1,2, ..,r  <  6  where  r  is  the  number  of  rigid  body  modes  possessed  by  the  FE  model.  The 
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T 

following  orthogonality  relations  are  assumed  to  hold  for  all  linear  (sub)structures:  $  M$  =  I 
and  =  diag|to2  j,  i.e.  the  modal  matrix  3>is  the  mass-normalized.  Also,  for  lightly 

damped  structures,  C  is  commonly  defined  to  be  proportional,  and  hence, 
C  =  M<I>  diag(2<;G))<I>TM  where  ^  is  the  dimensionless  damping  ratio  for  the  i*  mode.  This 

expression  is  valid  for  0  <  £  <  1  (underdamped).  For  an  underdamped  structure,  the  i*  damped 

natural  frequency  is  codi  . 


2.4  Inertial  and  Relative  Coordinate  Subsets 

Given  the  importance  of  base  excitation  to  this  work,  the  use  of  relative  coordinates 
instead  of  inertial  coordinates  will  allow  the  direct  application  of  prescribed  acceleration  to  the 
model,  eliminating  the  need  to  differentiate  a  base  displacement  or  velocity  time  history.  Of 
course,  if  inertial  displacements  or  velocities  are  to  be  recovered,  the  acceleration  time  history 
must  be  integrated. 

We  can  define  a  set  of  prescribed  ground  displacement  and  velocity  time  histories,  y^t), 
yi(t),  i=l,2,...,g.  We  also  define  g  subsets  xi  of  the  inertial  coordinates  x,  where  each  member  of 
each  subset  is  to  be  replaced  with  a  coordinate  relative  to  a  specific  ground  motion,  y^t).  These 
coordinate  subsets  are  denoted  Zj(t),  j  =  l,2,...,g,  where  each  relative  coordinate  set,  Zj  has  nj 
coordinates.  This  set  definition  allows  any  subset  of  inertial  coordinates  to  be  redefined  as 
relative  to  any  of  multiple  ground  motions  being  applied  to  the  model,  and  allows  the  retention 

g 

of  m  inertial  coordinates,  if  required.  Therefore,  the  length  of  the  z  vector  is  N  =  m  +  £  n ; . 

j=l 

The  coordinate  and  ground  motion  vectors  can  be  constructed  as  follows, 
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yi(t) 
yg(t). 

and  the  coordinate  transformation  between  the  exclusively  inertial  coordinate  vector  and  the 
mixed  inertial/relative  vector  is, 

x  =  z  +  Gy.  (2.3) 

Transforming  Eq.  (2.1)  provides, 

Mz  +  Cz  +  Kz  =  fe(t)  -  MGy  -  CGy  -  KGy.  (2.4) 

Since  CG=Cg  and  KG=Kg  (see  Eq.  (2.2)),  we  have, 

Mz  +  Cz  +  Kz  =  fir)(t)-MGy  =  fir)(t)-Mgy  (2.5) 

where  the  superscript  (r)  is  used  here  to  denote  that  elements  of  the  loading  vector  containing 
base  excitation  terms  have  been  “zero-ed,”  consistent  with  the  transformation  to  relative 
coordinates.  We  have  kept  the  base  acceleration  terms  separate  from  for  clarity,  although  we 
can  redefine  fe  =  f^  -  Mgy ,  recognizing  that 

Mz  +  Cz  +  Kz  =  fe  (2.6) 

is  essentially  the  same  equation  as  Eq.  (2.1),  but  where  the  use  of  the  symbol  z  emphasizes  that 
relative  coordinates  are  included. 
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2.5  Derivation  of  Governing  Equation  of  Nonlinear  Transient  Synthesis 

The  total  solution  for  (linear)  transient  response  can  be  written  in  terms  of  the 
convolution  integral, 

x(t)  =  xh(t)  +  JH(t-T)f(x)dT,  (2.7) 

0 

where  x  is  the  total  forced  response,  xh  is  the  homogeneous  solution,  f  is  the  excitation  vector, 
and  these  vectors  are  partitioned  according  to  the  previously  defined  sets  of  DOF,  e.g. 


x(t)  = 


xi(0 ' 

xc(t) 

xm(t) 

*b(0. 


f(t)  = 


fc(t) 

fm(t)  ' 

M'). 


(2.8) 


This  partitioning  is  implicit  in  all  matrix  equations  which  follow,  unless  otherwise  indicated.  The 
matrix  H  is  the  impulse  response  function  (IRF)  matrix,  any  element  of  which  can  be  written  as, 


Hij(t)=  £*frft+  | 

p=l  p=r+l  ®dp 


(^dp^)’ 


where  <t>f  is  the  i*  element  of  the  p*  eigenvector  of  the  substructure  prior  to  synthesis,  o)p  and  codp 
are  the  p*  undamped  and  damped  natural  frequencies,  respectively,  £p  is  the  p*  modal  damping 

ratio,  r  is  the  number  of  rigid  body  modes,  and  n<N  is  the  number  of  modes  required  for 
convergence.  The  number  of  elastic  modes  is  n-r.  Note  that  the  IRF  matrix  H  contains  elements 
from  all  substructures  involved  in  the  synthesis,  and  is  partitioned  as  described  above. 

We  can  decompose  each  excitation  component  into  an  externally  applied  portion  and  a 
component  due  to  synthesis,  as  follows.  The  i-set  DOF  are  by  definition  subject  only  to 
externally  applied  excitations,  so 

fi(t)  =  ff(t)  .  (2.10) 
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where  the  superscript  “e”  indicates  an  externally  applied  excitation.  The  c-set  DOF  are  subject  to 
externally  applied  excitations  as  well  as  coupling  reactions  hence, 

WO  =£(«)-*£(«)  (2.11) 

where  the  ~  overstrike  indicates  the  coupling  reaction,  R  is  a  Boolean  matrix  reflecting  the 
equilibrium  which  exists  between  the  coupled  DOF  [4],  and  the  *  superscript  indicates  a 
synthesized  (unknown)  quantity.  The  m-set  DOF  are  subject  to  externally  applied  excitations  as 
well  as  reactions  due  to  the  modifications  hence 

4>  W  =  fmW  -  fm(sm(t).*m(<))  =  (2.12) 

where  the  reactions  due  to  the  modifications  are  an  arbitrary  nonlinear  function  fm(t)  of  the 
synthesized  displacements  and  velocities  (accelerations  and  time).  The  b-set  DOF  are  subject  to 
externally  applied  excitations  as  well  as  excitations  due  to  prescribed  base  motion  y  acting 
through  an  arbitrarily  nonlinear  structural  element,  typically  involving  displacement-  and 
velocity-dependent  forces  only,  i.e. 

fb(t)  =  fb(t)-fb(xb(t)-y(t)’ib(t)-y(t))  =  fb(t)-fb(t)  (2.13) 

Therefore,  the  governing  equation  for  synthesis  is, 


x*(t)  =  x(t)- jH(t-x)f  *(t,T,x*(x),x*(T))dT  (2.14) 

0 

where  x(t)  contains  both  the  initial  displacement  and  response  due  to  externally  applied 
excitations, 

x(t)  =  x0  +  |H(t-T)fe(T)dT  (2.15) 

0 

and  f*(t)  are  the  synthesized  reactions  acting  on  all  DOF  sets. 
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(2.16) 


'  0  ' 

Rfc(t) 

Ot)  ‘ 

.  fb(t) . 

Equation  (2.14)  is  a  nonlinear  Volterra  integral  equation  of  the  second  kind,  and  is  the  central 
equation  of  this  work.  Direct  solution  is  possible  for  linear  problems;  for  nonlinear  problems 
iterative  solutions  are  required,  and  these  exploit  the  contractive  nature  of  the  integral  operators 
in  achieving  excellent  convergence  properties,  as  will  be  demonstrated. 


2.6  Solution  of  Governing  Equation  for  Synthesis: 

Uniqueness  and  Convergence 

In  order  to  investigate  issues  of  solution  uniqueness  and  convergence  of  Eq.  (2.14),  the 
method  of  Picard  Iteration  will  be  used  [9].  A  few  preliminary  concepts  and  definitions  from 
analysis  are  included,  without  proof,  for  completeness.  Proofs  may  be  found  in  [10,  Ch.  1]. 

We  define  the  following  metric  on  the  space  of  continuous  vector  functions  C(n)[0,t’]  as: 

=  max  ||x(t)  -  y(t)L  =  ||x  -  y||  (2. 17) 

t€[0,t  J 

where  x(t),  y(t)  e  C(n)[0,t’],  hence  defining  a  metric  space  C  =  (C(n),d).  A  (scalar)  sequence  (xn)  in 

the  metric  space  (C(n),d)  is  convergent  if  there  is  an  xe  C(n)  such  that  lim  d(xn,x)  =  0.  Then,  x 

n— »o° 

is  the  limit  of  (xn)  and  lim  xn  =  x ,  or  xn  — >  x .  A  sequence  (xj  in  the  metric  space  (C(n),d)  is 

n— 

Cauchy- convergent  if  for  every  8  >  0  there  is  an  N=N(e)  such  that  d(xm,xn)  <  s  for  every  m,n  > 

N.  The  metric  space  C(n)  is  complete  if  every  Cauchy  sequence  in  C(n)  converges,  i.e.  has  a  limit 
which  is  an  element  of  C(n).  This  distinct  and  important  characteristic  of  completeness  is  needed 
as  there  exist  spaces  in  which  Cauchy  sequences  do  not  converge.  Such  spaces  are  said  to  be 
incomplete.  It  can  be  noted  here  that  every  convergent  sequence  in  C(n)  is  a  Cauchy  sequence. 
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A  fixed  point  of  a  mapping  T:  C(n)  — »C(n)  of  a  set  C(n)  into  itself  is  an  x  e  C(n)  which  is 
mapped  onto  itself,  Tx  =  x.  Let  C  =  (C(n),d)  be  a  metric  space.  A  mapping  T:  C(n)  — »C(n)  is  a 
contraction  on  C®  if  there  is  a  positive  real  number  a  <  1  such  that  for  ajl  x,y  e  C(n), 

d(Tx,Ty)  <  ad(x,y)  (2.18) 

Central  to  the  effort  of  demonstrating  the  uniqueness  of  a  solution  to  Eq.  (2.14)  is  the 
Banach  Fixed  Point  Theorem.  Consider  again  the  non-empty  and  complete  metric  space  C  = 
(C(D),d)  and  let  T:  C(n)  — >  C(n)  be  a  contraction  on  C(n),  Then  T  has  precisely  one  fixed  point  [10, 

Sec.  5.1]. 

We  now  examine  the  nonlinear  Volterra  integral  equation, 

x*(t)  =  x(t)-jH(t-x)f(t,x,x*(x))dx  (2.19) 

o 

where  x(t)  is  bounded,  mt  <  ||x(t)||  .<  m2,  and  continuous  on  [0,t’]  as  it  is  the  solution  of  a  linear 
structural  dynamics  problem,  and  f(t,x,x*(x)j  is  assumed  to  be  continuous  (for  now),  and  hence 

bounded  on  the  rectangle  t  e  [0,t’],  x  e  [0,t’].  We  also  assume  that  ||F|t,x,x*(x)||<M  for  bounded 

x*(x),  i.e.  Xj<||x*(t)||<  xu.  We  note  that  x*(x)  can  be  reasonably  assumed  bounded  as  it  is  the 
(unknown)  solution  to  a  locally  nonlinear  structural  dynamics  problem  where  the  boundedness  of 
the  response  is  a  goal  of  design.  Note  also  that  f|t,x,x*(x)j  represents  the  action  of  some 

passive,  semi-active,  or  active  structural/mechanical  element,  and  hence  its  boundedness  is 
expected  due  to  physical  limitations. 

Following  [10],  consider  the  mapping  defined  by  the  successive  approximation  (Picard 
Iteration)  of  Eq.  (2.19), 

xn+i  (t)  =  x(t)  -  j H(t  - x)f  (t,  x,  X*  (x),  x*  (x))dx  a  t(x*  (x),  x*  (x))  (2.20) 

0 


24 


defined  on  the  space  of  continuous  functions,  C(n)[0,t’].  We  examine  the  following  elements  from 
a (bounded) sequence 

.••Ci(t),x;(t),x;+1(t)... 

along  with  the  metric,  Eq.  (2.17).  Our  first  goal  is  to  show  that  the  mapping  of  Eq.  (2.20)  is  a 
contraction,  and  to  this  end,  the  metric  is  used  to  measure  the  “distance”  between  inputs  xn_i(t) 
and  xn(t),  and  the  respective  outputs.  This  distance  is, 


d(xn+i(t)>Xn(t))  =  d(T(xJ(t),iJ(t)),T(xJ_1(tXii_1(t)))*... 


x(t)  -  f  H(t  -  x)f(t,x,x*  (x),<(x))dx  - 
0 


X(t)  -  }  H(t  -  T)f(t,T,Xn_1(x),X*_1(x))dT 
0 


(2.21) 


or, 


d(xn+l(0>Xn(t))  = 


L 

J  H(t  -  T){f(t,X,X^(x),X*(x))  -  f(t,X,Xn_1(x),Xn_1(x))}dX 


(2.22) 


or,  letting 


f  *(t,x)  =  f(t,x,x*  (x),x*  (x))  -  f(t,x,x*_1(x),x*_1(x)) , 


(2.23) 


yields  the  following  form  of  Eq.  (2.21) 


Xn+i(t)  -  x*  (t)||  =  d(x„+1(t),x^(t))  = 


jH(t-x)fn*(t,x)dx 


(2.24) 


In  order  to  be  able  to  relate  this  metric  of  the  outputs  to  the  corresponding  metric  of  the  inputs 
(and  demonstrate  the  contractive  property  of  the  mapping  Eq.  (2.20)),  the 


quantity 


'  djc 

xn(x)-xn_i(x)  must  be  isolated,  which  is  not  possible  in  the  form  of  Eq.  (2.24),  as 
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this  quantity  is  operated  on  by  the  function  f.  To  this  end,  we  impose  that  f  satisfy  the  Lipschitz 
condition, 

IfnCt^)!  <  L||x*  (x)  -  (2.25) 

for  |t,x,x„(x),  x„(x)j  and  (t,x,Xn_i(x),Xn_i(x)j  in  the  domain  of  f,  and  where  L  is  a  positive 
constant.  Note  that  the  contraction,  Eq.  (2.18)  is  a  Lipschitz  condition  with  L  =  a  <  1  [11]. 
Imposing  the  Lipschitz  condition  on  the  integral  of  Eq.  (2.24), 

|xn+l  (0  -  xn  (Ofl  -  k  j  ||H(t  -  t)||  •  Jx^  (t)  -  x„_j  (x)||dx  (2.26) 

0 

We  can  evaluate  Eq.  for  two  members  of  the  sequence,  Xj  (t)  and  X2(t) , 

11*3  (l) -  x2(t)||  -  ||H(t  -  x)||  •  Jx2(t)  -  x^  (x)||dx.  (2.27) 

0 


and,  ' 

|x5  (t)  -  X4  (t)||  <  L  J  ||H(t  —  x)||  -  ||x  J  -  xj  (x)||dx  (2.31) 

0 
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which  becomes, 


||x5(t)-X4(t)]<L3^-||H(t)||3|xu-x1|  (2.32) 

and  by  induction  the  following  can  be  established, 

|4+,(t)  -  4(t|  <  i^|l(Lt|H(t)||)n-1  (2.33) 

which  establishes  the  uniform  convergence  of  the  series, 

5  ||xn+l(t)-Xn(t|  (2.34) 

n=l  11 

with  no  restriction  on  L,  t,  or  ||H(t)||.  To  show  the  convergence  of  the  sequence  xn(t)  of  Eq. 

d|c  dfc 

(2.20)  to  the  solution  x  (t)  of  Eq.  (4.5),  we  write  xn(t)  using  a  telescoping  series, 

xnW  =  xl(t)  +  2  (xj+i(t)  -  xj (t)).  (2.35) 

Given  the  convergence  of  Eq.  (2.34),  taking  the  limit  of  Eq.  (2.35)  as  n  — >  °°  reveals  that 

lim  x„(t)  exists  for  all  1 6  [0,t'  ].  This  sequence  is  identical  to  the  sequence  generated  by  Eq. 
n— 

(2.20) .  We  assumed  that  f|t,T,Xn(t),Xn(t)jis  continuous,  and  hence  taking  the  limit  as  n  — » 
on  both  sides  of  Eq.  (2.20)  leads  to  the  conclusion  that  lim  xn(t)  =  x*  (t),  the  solution  to  Eq. 

n— 

(2.14). 

Although  we  have  shown  that  the  convergence  of  the  successive  approximation  is  not 
dependent  on  ||H(t)||,  it  is  worthwhile  to  evaluate  ||H(t)||  at  this  point.  Defining  the  p*  modal 
residue  matrix, 

r<p> =+<pyp>)T 

we  have  R  =  2  R®  =  M-1.  The  pUl  modal  impulse  response  function  is, 

P=1 
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-?v  .  /  \ 

e  sm((DdptJ 

^dp 


(2.36) 


and  hence  any  element  of  H,  Eq.  (2.9),  can  be  written  as, 


Hjj(t)  —  S  R[j  t+  £  Rg  hp(t) 

p=l  p=r+l 


(2.37) 


where  r  is  the  number  of  rigid  body  modes  and  n  is  the  total  number  of  modes.  An  upper  bound 
for  the  p*  modal  residue  is 

|r{-|  <  max^M-1),  p=l,2,..,n 
and  for  the  modal  impulse  response, 

max  |hp(t)|  <  — ^ — , 

0<t<t'  F  1  0)dp 


and  hence  the  norm  of  any  element  Hy  can  be  written  as 


f 


max  |Hij(t)|  < 

0<t<t’‘  J  1 


rf+  X  coi) 

p=r+l  p 


\ 

J 

y 

(2.38) 


Given  this  result,  we  find  that 


max  ||H(t)||  <  J  rt’+  Y  (Od1 
0<t<t'  p=r+l  p 


\ 


(2.39) 


where  J  is  the  number  of  columns  in  H. 
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3.  RECURSIVE  ALGORITHM  3: 

RECURSIVE  BLOCK-BY-BLOCK  CONVOLUTION 

The  algorithms  RA1  and  RA2  (Section  4)  are  similar  in  that  they  are  all  modal-based.  In 
these  algorithms,  the  recursion  is  obtained  from  the  use  of  a  modal  transition  matrix. 
Furthermore,  it  is  shown  that  these  algorithms,  which  are  inherently  unstable  in  their  explicit 
forms,  can  be  stabilized  by  reformulation  of  the  algorithms  in  implicit  forms,  specifically  as 
predictor-corrector  forms.  The  stability  achieved  through  the  use  of  a  single  corrector  step  is  only 
marginally  effective,  and  hence  we  develop  a  different  recursive  approach  for  the  solution  of  the 
governing  nonlinear  integral  equation,  Eq.  (2.14).  The  solution  method  to  be  developed  will 
exploit  the  exponential  convergence  result  of  Section  2.6  as  the  solution  method  is  iterative.  The 
algorithm  will  improve  greatly  upon  the  performance  of  the  basic  iterative  algorithm  reported  in 
[3],  by  the  development  of  a  recursive,  block-by-block  convolution  solution.  In  fact,  the  block- 
by-block  convolution  is  ideally  suited  for  the  calculation  of  structural  response  for  long  time 
records,  as  will  be  demonstrated. 

3. 1  Numerical  Quadrature  for  Nonlinear  Volterra  Integral  Equations 

The  numerical  solution  of  Eq.  (2.14)  typically  starts  with  a  discretization  of  the  equation  using 
some  quadrature  rule.  For  the  response  at  some  time  t  =  iAt  (t  =  0  =  OAt),  Eq.  (2.14)  becomes, 

x*(iAt)  =  x(iAt)-(At)“2wjH((i-j)At)f*(jAt)  (3.1) 

j=o 

where  we  have  abbreviated  the  general  nonlinear  force  as  f ,  a  and  (3  are  real  scalar  constants 
depending  on  the  quadrature  rule  chosen,  and  the  Wj  are  the  quadrature  weights.  For  example,  if 
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we  consider  the  simplest  of  quadrature  rules,  the  rectangular  rule  (for  a  purpose  to  be  made  clear 
below),  a  =  1,  |3  =  1,  and  w—l.  For  i  =  0,1,2,  Eq.  (3.1)  becomes , 

x*(0At)  =  x(0At)  (3.2) 

x’(lAt)  =  x(lAt)  -  At[H(lAt)f*(0At)  +  H(0At)f  *(lAt)]  (3.3) 

x'(2At)  =  x(2At)  -  At[H(2At)f*(0At)  +  H(lAt)f*(lAt)  +  H(0At)f*(2At)]  (3.4) 

and  we  note  that  H(t=0)  =  0,  yielding  the  correct  series  for  the  rectangular  rule.  These  equations 
can  be  rewritten  using  a  simplified  indexing  scheme,  i.e.  x(iAt)  =  x(i+l),  which  corresponds  to 

the  indexing  required  for  computer  programming,  e.g.  x(OAt)  =  x(l).  Using  this  indexing, 
Eqs.(3.2),  (3.3),  and  (3.4)  become,  for  i  =  1,2,3, 

x*(l)  =  x(l)  (3.5) 

x  *  (2)  =  x(2)  -  At[H(2)f  *  (1)  +  H(l)f  *  (2)]  (3.6) 

x  *  (3)  =  x(3)  -  At[H(3)f  *  (1)  +  H(2)f  *  (2)  +  H(l)f  *  (3)]  (3.7) 

where  again,  H(l)  =  0.  As  will  be  seen  below,  the  bracketed  terms  in  Eqs.(3.5),  (3.6),  and  (3.7) 
are  available  from  the  discrete  convolution. 

The  trapezoid  rule  and  Simpson’s  rule  are  more  commonly  used  quadratures  for  this 
application  [9,12].  The  performance  of  the  trapezoid  rule  in  the  solution  of  the  linear  synthesis 
problem  was  reported  by  Gordis  [2]. 
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3.2  Discrete  Convolution  and  Filter  Matrices 


We  define  the  basic  convolution  in  order  to  establish  a  notation  for  the  development  of 
the  block-by-block  convolution  which  follows.  The  convolution  of  two  vectors  x  and  y  is 
denoted  as  x*y.  The  discrete  convolution  of  x  and  y  is  given  by 


x*y  =  £x(n_k)y(k) 


(3.8) 


If  x  and  y  are  each  (n  x  1),  e.g. 


x  =  (x,  x, 

y =(y»  y2 


Xn-1  X„)T 

y».i  yn)T 


then  the  convolution  x*y  can  be  written  as  the  following  matrix- vector  product,  where  the 
matrix  is  Toeplitz,  constant  diagonal,  and  is  referred  to  as  a  filter  matrix  h(x)  [13]: 


z  =  x*y  =  h(x)-y  = 


>< 

o 

o 
_ 1 

r 

yi 

X,  Xj  0  .  : 

y2 

:  x2  X!  0  ••• 

• 

Xo-l  Xn-2  -  X1  0 

y«-i 

.  Xn  Xn-1  *„-2  •"  X2  Xl. 

.  yn . 

(3.9) 


For  example,  if  x  and  y  are  both  3  *  1,  we  have 


1 

>< 

o 

o 

_ 1 

yi 

x*y  = 

x2  Xj  0 

y2 

•  =  h(x)-y  =  < 

x2yi  +  x,y2  • 

x3  x2  Xu 

y3. 

.x3yi  +  x2y  2  +  xjj 

(3.10) 


where  the  filter  matrix  of  the  vector  x  is 


h(x)  = 


Xj  0  0 

x2  Xj  0 
x,  x,  x. 


(3-11) 
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and  the  elements  of  x  are  referred  to  as  filter  weights  [13].  Note  that  here  we  refer  to  h  as  an 
arbitrary  filter  matrix,  which  should  not  cause  confusion  with  the  use  of  the  symbol  H  to  refer  to 
the  impulse  response  function  (IRF)  matrix,  as  the  IRF  matrix  is  a  filter  matrix  as  well. 

From  a  comparison  of  Eqs.  (3.5),  (3.6),  (3.7)  with  Eq.  (3.10),  it  is  clear  that  the  discrete 
convolution  is  equivalent  to  the  use  of  the  rectangular  rule  for  numerically  integrating  a 
convolution  integral. 

We  now  define  a  delay  matrix  D  [13]  with  the  following  structure: 


. 

.  .' 

•  0 

0 

0 

0  • 

•  1 

0 

0 

0  • 

•  0 

1 

0 

0  • 

•  0 

0 

1 

0  • 

• 

• 

• 

• 

(3.12) 


where  the  dimension  of  D  is  consistent  with  the  length  of  the  vector  on  which  it  operates.  The 
matrix  D  produces  a  delay  in  time  by  one  sample.  For  example,  consider  the  3  by  1  vector  x. 


'0 

0 

o' 

V 

0‘ 

Dx  = 

1 

0 

0 

X, 

►  =  *< 

*1  ■ 

_0 

1 

0 

.x3. 

x2. 

(3.13) 


where  the  product  Dx  is  equivalent  to  the  vector  x  shifted  forward  in  time  (delayed)  by  one 
sample.  We  can  introduce  delays  of  an  arbitrary  number  of  samples  as  Dk.  The  product  Dkx 
produces  a  vector  equivalent  to  the  vector  x  but  delayed  by  k  samples. 

The  filter  matrix  h  is  equal  to  the  summation  of  powers  of  the  delay  matrix  multiplied  by 
the  filter  weights,  xs.  Alternatively,  the  columns  of  the  filter  matrix  h  are  each  products  of  powers 
of  the  delay  matrix  D  and  the  vector  x,  i.e.  the  k*  column  of  h  is  given  by  DS:.  The  filter  matrix 
of  a  vector  x  of  length  n  is  therefore. 
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(3.14) 


h(x)  =  Xxk ' °k  =[D°X  I>lx  •••  D"  2x  Dn"'x] 

k-0 


3.3  Block-by-Block  Convolution 

We  now  develop  the  block-by-block  (BBB)  convolution  of  two  vectors,  x  and  y,  i.e.  x*y. 
We  subdivide  the  entire  time  record  of  duration  T  seconds,  consisting  of  N  sample  points  (At  = 

T/N)  into  a  number  of  equally  sized  blocks,  or  subintervals,  i.e.  each  subinterval  contains  the 
same  number  of  sample  points.  We  will  subdivide  the  entire  record  into  “K”  blocks,  where  each 
block  consists  of  J  =  N/K  samples,  and  the  duration  of  each  block  is  JAt  seconds. 

It  is  important  to  emphasize  that  there  is  a  delay  of  J  samples  between  blocks.  For  the 
purpose  of  developing  the  BBB  algorithm,  we  will  need  to  extract  those  rows  of  a  vector 
corresponding  to  a  particular  block.  To  this  end,  we  define  the  following  row  extraction  matrix  r: 

"0  -010  - 

:  oio 

r  =  0  10  (3.15) 

:  0  10 
0  -  -01 

The  product  of  the  matrix  r  with  a  vector  x  is  the  subvector  of  x  consisting  of  the  rows  (samples) 
of  the  K*  (i.e.  last)  block,  i.e.  r  •  x  =  xK  where 

XK  =  [XJ(K-1)+1  "■  Xn-1  Xn] 

Using  the  delay  matrix  D,  we  can  define  a  matrix  which  extracts  the  rows  of  the  k*  block,  where 
k=  1,2,...,K. 

'0  -  0  10--  •••  0" 
i  0  10  i 

rk  =  r  •  Dk  =  0  1  0  (3.16) 

:  0  10  ; 

0  -  -"010  •••  0 
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The  matrix  equation,  Eq.  (3.9),  can  be  written  in  a  block-partitioned  form  as  follows.  We  can 
write  the  k*  subvector  of  z,  i.e.  zk,  as, 

zk=Srk'h(x)^-ym  (3.17) 

m=l 

where  ym  =  rm  •  y .  It  is  important  to  note  that  the  block  filter  matrix  rk  •  h(x)  •  rj  need  never  be 
formed,  as  the  following  relations  hold: 

T  jx(l:J)  if  k  =  m 

rk  Mx)  rm  |x|((k-m  +  l)J-(2J-2)):(k-m  +  l)j)  ifk>m  (-•  '-) 

where  x(p:q)  indicates  the  subvector  of  x  consisting  of  elements  p  through  q. 

3.4  Performance  Comparison  -  Standard  and  Block-by-Block  Convolution 

A  traditional  (single-block)  convolution,  for  sufficiently  long  records  of  length  n,  is  most 
efficiently  computed  using  the  FFT,  yielding  a  total  number  of  floating  point  operations 
(FLOPS)  proportional  to  n*log2(n).  The  computing  language,  MATLAB  provides  a  built-in 
function  for  convolution  which  uses  FIR  filters  for  the  calculation,  and  yields  total  FLOPS 
proportional  to  n2.  As  we  are  here  interested  in  comparing  the  performance  of  the  BBB  algorithm 
with  the  traditional  single-block  convolution,  the  use  of  the  MATLAB  function  will  provide 
much  convenience  with  no  loss  in  the  ability  to  compare  algorithms.  Of  course,  if  one  were  to 
optimize  a  convolution  algorithm  for  large  n,  the  FFT  approach  would  be  best. 

Assuming  one  iteration  for  each  diagonal  block,  the  number  of  FLOPS  for  the  BBB 
algorithm  is  given  by: 

FLOPS  =  K(2J2-J)  +  -(K2-K)(4J2-4J  +  1)  (3.19) 

2 

which  yields  an  optimum  number  of  blocks, 

K„=|+2N  (3.20) 


34 


which  yields  a  solution  which  is  clearly  not  realizable.  The  optimal  number  of  blocks  calculated 
is  greater  than  the  total  number  of  samples  N,  and  is  a  non-integer  number  of  blocks.  What  is 
useful  about  this  solution  is  that  is  indicates  that  the  FLOPS  required  by  a  BBB  convolution 
decreases  monotonically  with  increasing  block  number.  This  is  shown  in  Figure  2,  which 
compares  the  FLOPS  required  by  a  standard  convolution  to  the  BBB  convolution  for  varying 

total  number  of  samples,  N,  and  for  different  numbers  of  blocks. 
xIO8 

12 

10 

8 

FLOPS  6 

.  4 

2 

0 

Figure  2.  FLOPS  required  versus  convolution  vector  for  varying  number  of  blocks 

However,  if  we  compare  actual  compute  time  (using  MATLAB),  we  see  that  there  is  a  point  at 
which  increasing  the  number  of  blocks  results  in  increased  compute  times,  as  the  computing 
“overhead"  associated  with  increased  block  numbers  outweighs  the  decrease  in  computing  time 
due  to  the  reduction  in  FLOPS  required.  This  is  shown  in  Figure  3. 


r111 . i . i . ) . » 
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Figure  3.Computer  time  (sec)  required  versus  convolution  size  for  varying  number  of  blocks 


3.5  Recursive  Block-by-Block  Iteration  Solution 

Before  discussing  the  recursive  block-by-block  iteration,  we  outline  the  basic  iteration 
algorithm  for  the  solution  of  Eq.  (2.14).  In  the  algorithm  which  follows,  it  is  implied  that  the 
vector  x  is  partitioned  consistently  with  Eq.  (2.16),  with  the  alteration  that  the  partition 
associated  with  the  “i”  coordinates  has  been  deleted.  Only  those  coordinates  x*  directly  involved 
in  the  synthesis,  i.e.  those  coordinates  subjected  to  forces  of  synthesis,  are  included  in  the 
iteration.  The  “i”  set  coordinate  responses  are  calculated  by  a  direct  convolution  of  the  associated 
IRF  with  the  (converged)  forces  of  synthesis,  which  result  from  the  iteration.  The  coordinate  set 
involved  in  the  synthesis  is  the  defined  by  the  set  union 

s  =  mucub  (3.21) 

where  s  denotes  the  synthesis  set.  The  IRF  matrix  is  therefore  more  fully  denoted  as  H^t). 
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For  clarity  of  presentation,  the  time  dependence  and  asterisk  (•)*  indicating  a  synthesized 
quantity  will  be  dropped. 

Basic  Iteration : 

•  Initialize: 

j«=l 

fsj<=l 

•  While  xj+1  *  x] 

•  xj+1  <=xs-Hss*f(x],xj,y) 

•  j<=j  +  l 

•  Converged  forces  of  synthesis: 

C  <=f(x',x>,y) 

•  Solution  for  i-set  responses: 
x*=Xi-His*f; 

We  will  now  expand  this  algorithm  to  incorporate  the  recursive,  block-by-block 
approach.  The  algorithm  is  recursive  in  that  the  iteration  performed  for  block  “k”  makes  use 
of  the  already  converged  forces  of  synthesis  f'  for  prior-time  blocks  k-1,  k-2,  etc,  where  for 
the  sake  of  clarity,  the  “s”  subscipt  has  been  dropped.  As  will  be  described,  only  those  forces 
of  synthesis  at  the  current  block  are  included  in  the  iteration,  as  prior  block  synthesis  forces 
are  converged.  We  will  denote  the  responses  and  forces  for  the  k*  block,  and  at  the  jth 
iteration,  as  xJkand  fk.  The  IRF  filter  matrix  relating  the  k*  response  block  and  the  m*  input 
block  is  denoted  as  H^,,  and  is  given  by 

Hkm  =  •  H  •  r*  (3.22) 

There  are  K  blocks,  k  =  1,2,...,K,  and  each  block  is  of  length  J  (samples).  We  will  make 
use  of  Eq.  (3.22)  to  symbolically  denote  the  IRF  matrix  blocks,  while  keeping  in  mind  that  in 
practice  these  matrices  need  never  be  formed.  What  is  formed  in  practice  are  the  partitioned 
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vectors  from  which  these  IRF  blocks  are  constructed,  as  given  by  Eq.  (2.1).  The  iteration  for 
the  k*  block  is  given  by. 


(3.23) 

m=t 

We  can  now  construct  the  algorithm. 

Recursive  Block-by-Block  Convolution  Algorithm 

•  Initialize: 
j^=l 

fj  «=  1  (over  all  blocks) 

•  Do  k  =  1 :  K 

•  While  xJk+1  *  xJk 

m=l 

•  ff1 

•  j<=j+l 

•  End  While 

•  Converged  forces  of  synthesis: 

^k  ^ 

*, 

•  End  Do 

•  Solution  for  i-set  responses: 

x*=Xi-His*fs* 
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3.6  Example  of  Block-by-Block  Synthesis 

The  following  example  is  taken  from  [3],  in  order  to  demonstrate  the  improvement  in 
performance  of  the  current  block-by-block  algorithm  relative  to  the  basic  iteration  algorithm 
reported  in  [3]. 

This  example  demonstrates  the  use  of  the  nonlinear  synthesis  procedure  in  a  nonlinear 
base  isolation  problem  subjected  to  a  prescribed  transient  base  displacement  excitation.  An 
idealized  “deck”  model  is  isolated  by  four  nonlinear  springs/dampers  located  at  the  four  comer 
nodes,  as  shown  in  Figure  4.  A  piece  of  equipment  is  mounted  elastically  on  the  deck,  modeled 
using  a  single  lumped  mass  and  linear  spring.  The  deck  is  modeled  using  four-noded 
quadrilateral  elements.  The  deck  is  a  square  steel  plate,  10  inch  on  a  side,  and  the  plate  thickness 
is  0.125  inch.  The  lumped  mass  is  taken  as  5.5%  of  the  total  plate  mass.  The  linear  spring 
stiffness  is  1000  lbf/in. 


Figure  4.  Isolated  deck  with  equipment  (lumped  mass)  mounted  with  linear  spring 
The  excitation  used  is  a  prescribed  “blast”  base  motion,  shown  in  Figure  5. 
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Figure  5.  Transient  base  motion 

We  now  present  the  results  obtained  from  the  current  recursive  block-by-block  algorithm. 
In  this  example,  all  modes  up  to  12KHz  were  calculated  (99  modes).  The  synthesis  was 
performed  using  all  modes  up  to  4KHz.  The  isolators  used  in  this  example  are  described  by  the 
following  equation: 

f(x-y,x-y)  =  k(x-y)  +  k3(x-y)3  +  c2(x-y)2  (3.24) 

where  k  is  the  linear  stiffness  coefficient,  k3  is  the  cubic  stiffness  coefficient,  and  c2  is  the 
quadratic  damping  coefficient.  In  the  example  problem,  the  isolator  parameters  have  the 
following  values: 


k  =  100  lbf/in  k3  =  201bf/in3  C2  =  0.001  lbf-sec/in2 

The  synthesis  was  performed  using  the  following  parameters: 


Sample  length: 

5.0e-5  sec 

Modes  Retained: 

38 

Number  of  Subintervals: 

8 

Max  mode  frequency: 

3,702.9  Hz 

Samples/subinterval : 

101 

Modal  Damping: 

0.0 

End  time: 

0.04  sec 
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The  Direct  FE  solution  used  the  same  sample  length  and  end  time  as  the  synthesis. 

In  Figure  6,  the  vertical  displacement  response  of  a  comer  node  of  the  plate  model  is 
shown,  along  with  the  corresponding  response  quantity  as  calculated  using  the  commercial 
Direct  FE  solution.  Keeping  in  mind  the  relative  times  required,  7  min  54  sec  for  the  synthesis 
(including  the  solution  for  99  modes  up  to  12KHz),  versus  30  min  15  sec  for  the  Direct  FE 
solution,  its  clear  that  the  synthesis  provides  a  very  accurate  solution. 

The  corresponding  comparison  of  the  velocity  for  the  same  node  is  shown  in  Figure  7.  While 
there  is  some  disparity  in  the  two  solutions  for  velocity  at  early  times,  these  differences  do  not 
compromise  the  solution  after  the  early  times.  This  is  due  to  the  fact  that  the  convergence  of  the 
numerical  method  is  relatively  insensitive  to  errors  in  the  starting  values;  the  effect  of  the  starting 
errors  on  convergence  is  attenuated  by  the  factor  At  [9]. 

In  Figure  8,  the  vertical  displacement  response  of  the  lumped  mass  is  compared.  It  should 
be  noted  that  the  synthesis  solution  for  the  response  of  the  lumped  mass  is  a  straight  convolution 
using  the  nonlinear  forces  which  are  the  direct  result  of  the  synthesis.  The  isolated  mass  is  a 
model  DOF  not  subjected  to  forces  of  synthesis  (i.e.  nonlinear  isolator  forces),  and  hence  this 
DOF  is  a  member  of  the  coordinate  set  “i”  rather  than  the  “c”  coordinate  set. 
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Figure  7.  Deck  comer  vertical  velocity  response 


3.6. 1  Comparison  of  Computing  Times  Required:  Synthesis  vs.  Direct  FE 

In  the  tables  which  follow,  a  comparison  is  made  of  the  computing  times  required  for  a 
complete  transient  analysis  as  performed  using  the  synthesis  and  using  commercial  nonlinear 
direct  transient  analysis.  The  direct  transient  analysis  does  not  use  the  modes,  and  hence  the  total 
time  for  the  direct  FE  solution  does  not  include  a  modal  solution.  The  synthesis  does,  however, 
require  a  modes  database  for  each  substructure,  and  this  solution  is  performed  using  the 
commercial  FE  program.  The  time  required  for  a  single  complete  transient  analysis  using 
synthesis  therefore  includes  the  time  for  the  modes  solution  and  for  the  synthesis  itself.  As  is 
seen  from  the  tables,  all  subsequent  analyses  require  only  the  synthesis  itself,  and  hence  an 
enormous  savings  in  time  is  realized,  as  compared  with  the  direct  FE  solution,  which  must  be 
repeated  in  its  entirety  for  each  analysis.  All  of  the  calculation  times  reported  in  each  table  were 
generated  using  a  common  computer  which  ran  both  the  commercial  FE  program  and  the 
previously  reported  synthesis  algorithm,  written  in  MATLAB. 
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Table  1.  Summary  of  Solution  Times  Required:  Recursive  Block-by-BIock  Algorithm 


Computing  Time  Required 


Calculation  Performed 

Using  Synthesis 

Using  Direct  FE 

Normal  Modes: 

(All  modes  to  12KHz) 

7  min  47  sec 

not  required 

Transient  Response: 

0  min  7  sec 

30  min  15  sec 

Total  Time  Required  - 
Single  Analysis 

7  min  54  sec 

30  min  15  sec 

Total  Time  Required  - 
Subsequent  Analyses 

0  min  7  sec 

30  min  15  sec 

The  following  table.  Table  2,  taken  from  Reference  [3],  shows  the  comparison  of  the 
previously  reported  algorithm  described  in  [3]  with  a  commercial  direct  nonlinear  transient  FE 
solution.  This  table  is  included  here  to  demonstrate  the  improved  performance  of  the  new 
recursive  block-by-block  algorithm  as  compared  with  the  algorithm  of  Reference  [3].  As  can  be 
seen  from  the  times  reported  in  Table  2,  for  a  single  analysis,  the  original  algorithm  required 
approximately  60%  less  computer  time  as  compared  with  the  commercial  nonlinear  direct 
transient  analysis.  The  new  algorithm  requires  approximately  74%  less  computer  time  than  the 
commercial  nonlinear  direct  transient  analysis.  However,  the  improvement  in  computing  time  for 
the  synthesis  alone  is  dramatic;  the  new  algorithm  requires  approximately  7  seconds  to  perform 
each  subsequent  nonlinear  transient  analysis.  This  fast  reanalysis  capability  will  facilitate  the  use 
of  numerical  optimization  techniques  for  nonlinear  isolation  design. 
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Table  2.  Summary  of  Solution  Times  Required:  Previously  Reported  Algorithm  Ref.  [3] 


Computing  Time  Required 

Calculation  Performed 

Using  Synthesis 

Using  Direct  FE 

Normal  Modes: 

(All  modes  to  12KHz) 

13  min  24  sec 

not  required 

Transient  Response: 

4  min  15  sec 

43  min  41  sec 

Total  Time  Required  -Single 
Analysis 

17  min  39  sec 

< 

43  min  41  sec 

Total  Time  Required  - 
Subsequent  Analyses 

4  min  15  sec 

43  min  41  sec 

3.7  Summary  of  Recursive  Algorithm  3:  Block-by-Block  Convolution 

The  block-by-block  convolution  algorithm  provides  a  dramatic  decrease  in  computer  time 
required  for  convolutions.  As  it  was  shown  that  the  discrete  convolution  is  identical  to  a 
numerical  integration  using  the  rectangular  rule,  the  block-by-block  algorithm  provides  the  same 
decrease  in  computer  time  required  for  the  solution  of  the  governing  nonlinear  integral  equation. 
While  the  time  savings  for  a  single  analysis,  relative  to  a  direct  nonlinear  transient  analysis  is 
large  (approximately  74%  for  the  example  presented),  the  time  savings  for  each  subsequent 
analysis  is  extremely  large,  as  is  summarized  in  Table  1.  This  extremely  fast  reanalysis  will 
facilitate  the  use  of  numerical  optimization  for  locally  nonlinear  structures. 

Currently,  each  block  is  integrated  using  the  convolution  (i.e.  rectangular  rule).  However, 
it  is  a  straightforward  matter  to  implement  the  Trapezoid  rule,  which  will  provide  an  increase  in 
accuracy.  The  algorithm  lends  itself  to  parallel  computation,  and  to  variable  step  size/block  size 
as  well.  In  its  simplest  form,  the  step  size  for  each  block  can  be  established  independently  of  the 
other  blocks.  For  example,  a  short  initial  block  with  a  small  step  size  can  be  used  to  minimize 
starting  errors. 
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4.  RECURSIVE  TRANSITION  MA  TRIX  ALGORITHMS 

The  algorithms  developed  in  this  section.  Recursive  Algorithms  1  and  2  (RA1  and  RA2), 
were  developed  prior  to  Recursive  Algorithm  3  (RA3),  which  is  presented  in  Section  3  of  this 
report.  Algorithms  RA1  and  RA2  are  based  on  modal  transition  matrices,  and  hence  are 
fundamentally  different  form  RA3.  As  was  discussed  in  Section  1.2,  these  algorithms  did  not 
meet  the  goal  stated  in  Section  1.1,  but  are  documented  for  completeness,  and  for  the  potential 
benefit  of  the  developments  in  other  applications,  specifically  with  regard  to  RA2.  We  begin 
with  some  background. 

Many  approaches  to  the  transient  analysis  of  locally  nonlinear  problems  are  direct,  in  that 
they  are  based  on  the  physical-coordinate  (not  modally  transformed)  second-order  differential 
equations  [14,  Sec.  2.C.1], 

Mx  +  fnl(x,x)  =  fe  (4.1) 

where  x  is  a  vector  of  “N”  nodal  DOF,  and  the  ('),(")  indicate  the  first  and  second  time 

derivatives.  M  is  an  N  by  N  symmetric  positive-definite  mass  matrix,  feis  a  vector  of  externally 
applied  loads,  and  ^represents  the  loading  which  is  a  nonlinear  function  of  the  displacement  and 
velocity  vectors.  In  the  case  where  the  damping  and  elastic  forces  have  linear  portions  separable 
from  the  nonlinear  portions,  Eq.  (4.1)  can  be  written  as 

Mx  +  Cx  +  Kx  =  fe  +  fnl  (x,  x)  (4.2) 

where  C  is  a  symmetric  positive-semidefinite  damping  matrix,  and  K  is  a  symmetric  positive- 
semidefinite  stiffness  matrix.  These  three  matrices,  along  with  M,  represent  the  linear  portion  of 
the  model.  In  Eq.  (4.2),  the  vector  fn|(x,x)  represents  the  forces  imposed  on  the  linear  portion  of 
the  model  by  nonlinear  elements  in  the  model.  It  should  again  be  noted  that  we  are  limiting  our 
attention  to  localized  nonlinear  components  which  have  no  internal  DOF.  Such  a  component 
could  not  be  represented  by  Eq.  (4.2),  but  rather  by  the  following  system  of  equations, 
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Mx  +  Cx  +  Kx  =  fe  +  fnj  (x,x,z,z) 
z  =  gni(x,x,z,z) 


(4.3) 


(4.4) 

where  z  is  an  n*l  vector  of  DOF  internal  to  the  nonlinear  components,  and  g^is  the  nonlinear 
vector  function  representing  the  nonlinear  components.  Note  that  we  have  excluded  nonlinear 
inertia  terms  for  reasons  of  clarity,  not  necessity.  As  discussed  in  [14,  Sec.  2.C.1],  Eqs.  (4.1) 
through  (4.4)  can  describe  spatially-discrete  models  with  nonlinear  elastic  and/or  visco-elastic 
elements,  as  well  as  geometrically  nonlinear  elements.  However,  phenomena  such  as  plasticity 
and  visco-plasticity  are  not  addressed  by  these  equations  due  to  the  history-dependence  of  the 
associated  internal  forces  of  these  phenomena.  A  large  variety  of  algorithms  exist  for  the  solution 
of  Eqs.  (4.1)  through  (4.4)  and  are  well-documented.  See  Hughes  [15]  for  a  detailed  summary  of 
many  current  methods.  A  commonly  used  algorithm  in  structural  dynamics  is  the  average- 
acceleration,  unconditionally  stable  Newmark  method  [16]. 

It  is  also  possible  to  address  these  problems  in  first-order  form,  i.e.  . 

y  =  Ay  + Bfnl(y,t)  (4.5) 

where  y  is  a  (2N  x  1)  vector  contain  both  displacements  and  velocities,  A  is  the  (2N  x  2N)  non- 
symmetric,  positive  semi-definite  system  matrix,  B  is  the  (2N  x  b)  input  matrix  (for  “b”  inputs), 
and  f„i(y,t)  is  the  (b  x  1)  vector  of  nonlinear  loads  associated  with  this  equation.  More  will  be  said 
later  about  the  first-order  form  for  locally-nonlinear  transient  structural  dynamics. 

A  disadvantage  of  direct  methods  is  the  requirement  to  retain  all  DOF  of  the  model, 
regardless  of  whether  the  response  of  these  DOF  are  of  interest.  In  contrast  to  the  direct  methods, 
the  use  of  modal  methods  in  the  solution  of  the  differential  equations  (4.2)  and  (4.5)  can 
eliminate  the  need  to  retain  all  model  DOF,  but  a  disadvantage  with  such  a  modal  approach  is 
that,  strictly  speaking,  the  number  of  retained  modes  required  can  be  assessed  only  after  (two) 
solutions  have  been  performed,  thereby  demonstrating  modal  convergence. 
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The  direction  of  this  work  is  the  development  of  a  solution  method  of  Eq.  (4.12)  which  is 
designed  to  have  an  implicit  and  exact  model  reduction,  in  that  only  those  model  DOF  directly 
.  associated  with  externally  applied  and  nonlinear  forces  need  be  retained  in  the  analysis,  and  to  be 
numerically  robust  with  excellent  convergence  properties.  This  integral  equation  formulation  for 
transient  structural  synthesis  [2,3],  has  been  shown  to  provide  significant  reductions  in  computer 
time  as  compared  with  direct  integration  using  a  well-known  commercial  finite  element  program 
[3],  due  to  this  implicit  reduction. 

We  will  describe  two  transition  matrix  algorithms  developed.  We  discuss  the  motivation 
for  these  approaches,  why  these  algorithms  are  not  suitable  for  the  purpose  at  hand,  and  what 
insights  and  results  of  value  were  gained  by  the  development. 

4. 1  Recursive  Solutions  of  First-Order  Differential  Equations 

We  begin  with  a  brief  review  of  recursive  convolution  solutions  to  first-order  differential 
equations,  due  to  the  development  of  algorithms  RA1,  RA1PC,  RA2,  and  RA2PC.  The  following 
brief  review  can  be  found  in  many  texts  in  more  detail;  see  for  example,  Meirovitch  [17]. 

A  recursive  solution  is  attractive  due  to  its  reduced  storage  requirements,  which  are  non¬ 
trivial  when  processing  large  FE  models,  and  due  to  the  reduction  in  compute  time  which  can  be 
obtained,  depending  both  on  the  algorithm,  and  the  structure  representation  employed. 

The  total  solution  of  the  linear  first-order  equation, 

y  =  Ay  +  Bf(t)  (4.6) 

is  given  by, 

x(t)  =  eAtx(0)  +  JeA(t_x)Bf(T)dT  (4.7) 

0 

Evaluating  the  solution,  Eq.  (4.7),  at  sample  points  k  and  k+1,  gives  respectively 
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(4.8) 


kT 

x(kT)  =  eAkTx(0)  +  J  eA(kT'T)Bf(x)dx 
0 

kT+T 

x(kT  +  T)  =  eA^kT+T^x(0)  +  J  eA^kT+T_T^Bf(x)dx  (4.9) 

0 

A  fundamental  assumption  required  to  produce  an  explicit  recursive  (discrete)  form  of  Eq.  (4.7) 
is  that  the  force  f  is  constant  over  the  sample  interval,  i.e.  f(x)=fkT)  for  kT  <  x  <  (k+l)T.  With  this 

A  X 

assumption,  and  making  the  change  of  variable,  cr  =  kT  +  T-x,  and  defining  <E>  =  e  and 
T 

r  =  J  eAodcB  =  A-1(eAT  -  IjB,  the  following  discrete  recursive  form  of  Eq.  (4.7)  is  found: 

0  '  7 

x(k+l)  =  ^x(k)+rf(k)  (4.10) 

The  advantage  in  the  use  of  Eq.  (4.10),  as  opposed  to  Eq.  (4.7),  is  as  follows.  Considering  a 
single-input,  single-output  response  calculation  over  K  sample  points,  the  convolution  integral  in 
Eq.  (4.7)  requires  on  the  order  of  k2  flops,  and  at  each  sample  point  k,  all  prior  values  of  f  and 

are  required,  k-1,  k-2,  k-3,  etc.  The  recursive  form,  i.e.  Eq.  (4.10)  requires  only  the  data  at 
sample  point  k-1  to  calculate  the  response  at  sample  point  k.  This  constitutes  a  large  saving  in 
both  compute  time  and  data  storage  required.  While  these  advantages  are  significant,  the 
calculation  of  the  matrix  exponential  <J>  is  required.  As  the  system  matrix  A  is  2n  by  2n,  the 

explicit  computation  of  is  not  practical  for  large  models.  Modal  calculations  of  <I>  are  possible, 

but  require  the  calculation  of  the  left  eigenvectors  of  A.  This  is  not  an  attractive  option  either, 
given  that  most  commercial  finite  element  programs  do  not  provide  for  this  calculation,  although 
it  is  possible  in  NASTRAN  using  the  DMAP  language.  Our  motivation  is  therefore  to  develop  a 
recursive,  transition-matrix  solution  algorithm  solution  which  has  an  inherent  implicit  exact 
model  reduction,  and  the  above  described  advantages  of  a  recursive  procedure.  We  also  stipulate 
that  any  such  method  not  require  the  calculation  of  exp(AT),  the  matrix  exponential  of  the  2N  x 
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2N  system  matrix.  Our  intermediate  goal  here  is  the  development  of  a  recursive,  transition  matrix 
approach  based  on  the  second-order  differential  equation  model,  i.e.  Eq.  (2.1). 

A  recursive  predictor-corrector  algorithm  for  locally  nonlinear  earthquake  isolation  was 
developed  by  Inaudi  and  De  La  Llera,  and  implemented  in  the  IN  ADEL  program  [7].  This 
algorithm  represents  the  structure  with  a  state-space  model,  i.e.  Eq.  (4.5).  The  algorithm 
developed  in  [7]  is  rendered  implicit  by  the  use  of  non-constant  approximations  for  the  force  f, 
which  is  defined  quite  generally.  While  this  work  contributes  to  the  motivation  for  exploring  a 
recursive  transition  matrix  approach,  as  stated  above,  a  goal  of  the  current  work  is  to  avoid  the 
calculation  of  the  transition  matrix  (exponential  of  the  system  matrix),  due  to  the  attendant 
computational  demand. 

4.2  Recursive  Algorithm  1:  First-Order  Modal  Transition  Matrix 

The  modal  transformation,  x=$>q,  of  Eq.  (2.1),  for  proportional  damping,  produces  a  set 
of  uncoupled  modal  differential  equations  of  the  form, 

q  +  diag(2qiO)i)q  +  diag(to?)q  =  f  (4.11) 

where  the  modal  force  is  f(t)  =  4>Tf(t),  containing  mass-normalized  mode  shapes,  and  where 
(~)  indicates  a  modal  quantity.  In  the  developments  which  follow,  the  force  f(t)  may  be 
interpreted  to  be, 

f(t)  =  f(t,x  *  (t),x  *  (t)) 

The  set  of  (at  most)  N  equations,  Eq.  (4.11),  is  comprised  of  two  basic  forms.  The  rigid  body 
modes  (tDj  =  0)  are  described  by 

qi=fi  i=l,2,..,r  (4.12) 
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where  r  is  the  number  of  rigid  body  modes  (RBM)  possessed  by  the  structure,  and  excluding 
mechanisms,  r  <  6.  The  solution  to  Eq.  (4.12)  is  given  by, 

t  x  t 

qi(0  =  J I fi(°)dGdT  +  qio  +  qiot  =  J(t  - T)fj(x)dT  +  q=  +  q:  t  (4.13) 

00  0 

The  elastic  modes  (co5  >  0)  are  described  by: 

qj +2^0)^ +0)^=^  (r+1)  <  i  <  N  (4.14) 

where  the  number  of  elastic  modes  retained  is  typically  much  less  than  N,  the  total  number  of 
nodal  DOF. 

We  may  also  write  the  solution  to  Eqs.  (4.12)  as  a  convolution  integral, 

t  ^ 

qi(t)  =  Jhi(t-T)fi(T)dT  +  qio.+  qiot,  (4.15) 

0 

where  qi()  and  qi()  are  the  modal  initial  conditions,  and  where  hj(t)  is  the  modal  impulse 
response  function  (IRF).  As  seen  from  Eq.  (4.13),  the  IRF  for  an  RBM  is, 

hi(t)  =  t,  i=l,2,..,r,  (4.16) 

and  for  an  elastic  mode,  the  modal  IRF  is  given  by, 

hi(t)  =  — e“?iCOitsin(o)dit),  (r+l)<i<N.  (4.17) 

“di 


We  will  construct  the  modal  transition  matrices  for  the  rigid  body  and  elastic  modes,  and 
the  associated  recursive  algorithm,  based  on  the  physical  (nodal)  and  modal  state  vectors  defined 
as  follows: 


51 


Based  on  Eqs.  (4.15),  (4.16),  and  (4.18),  the  solution  to  the  rigid  body  mode  equation,  Eq.  (4.12), 
can  be  written  as 

4i(t)  =  Vrb(t)qio  +/Yrb(t-T)Bfi(T)<fc  (4.19) 

0 

where  B  =  (0,1  )T  and  where  the  matrix  ^(t)  is  defined  as, 

^b(t)=  ‘  •  (4.20) 

In  order  for  a  matrix  Y  to  be  a  transition  matrix,  it  must  satisfy  the  following  property: 

V(t  +  At)  =  Y(t)Y(At)  =  Y(At)Y(t)  (4.21) 

It  is  easily  shown  that  satisfies  this  property,  and  hence  *Frb  is  the  RBM  transition  matrix,  as 
will  be  shown  here.  We  can  write  Eq.  (4.19)  for  a  time  one  sample  later,  i.e.  at  t+At, 

t+At 

qj(t  +  At)  =  Yrb(t  +  At)qio  +  J  Yrb(t  +  At-xjBfj^dT,  (4.22) 

0 

and  as  was  done  to  obtain  Eq.  (4.10),  in  conjunction  with  the  fact  that  Yrb  satisfies  Eq.  (4.21),  we 
find  the  recursive  form  of  Eq.  (4. 19), 

qi(t  +  At)  =  Trb(At)qi(t)  +  rrb(At)fi(t)  (4.23) 

where 

At  r i  ,5Ati 

rrb(At)  =  J  Yrb(a)daB  =  At  B.  (4.24) 

o  Lu  1  . 

If  we  define  the  k*  time  step  tk  =  kAt,  denote  a  quantity  evaluated  at  the  km  step  as  «(k),  and 
recognize  that  and  rrb  are  constants  for  a  given  sample  length  At,  Eq.  (4.23)  can  be  written  as, 

qi((k  +  l)T)  =  'Frb(At)qi(kT)  +  rrt(At)fi(kT)  i=1.2 . r  (4.25) 
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Focusing  on  an  elastic  mode,  the  state-space  equation  of  motion  for  an  elastic  mode  is, 


qi=Aiqi+Bfi 


(4.26) 


and  the  recursive  form  of  the  solution  to  Eq.  (4.26)  is. 


q{k+1)  =  'Fiq[k)  +  rA  (r+1)  <  i  <  N. 


(4.27) 


and  where  ^  =  eA‘ At ,  I\  =  Aj  -  l)B,  B=[0  1]T,  and  the  mode  system  matrix  is, 


-1 


Ai  = 


0 


1 


L-®i  -2Si®iJ 


(4.28) 


The  explicit  form  of  RA1,  based  on  Eq.  (4.25)  is  given  here.  The  number  of  modes  retained  is 
indicated  by  “p”  where  typically,  p«N. 


Recursive  Algorithm  #1  -  Explicit  Form 

Evaluate  Modal  Transition  Matrices  W  and  Input  Matrices  T 

k  :=  0 
for  i  =  1 :  p 

fp>;=*Tf(kT,x<k>',K«') 


if  i  <  r  then 

else 

end 


a<k+1)  = 


=>rrbq!k)+rrt,fjk» 
'FiqjK'+rifI(k) 


q(k+l)  _  vp.^(k) 


end 

$(k+l)  =  ^>q(k+1) 
k  :=  k  +  1 
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4.3  Recursive  Algorithm  2: 2nd-Order  Complex  Modal  Transition  Matrix 

This  algorithm  is  based  on  a  newly  developed  transition  matrix  constructed  from  the 
solution  to  the  second-order  elastic  mode  differential  equation.  The  key  to  the  development  of 
such  a  transition  matrix  is  to  write  the  solution  of  the  second-order  equations  with  a  leading 
matrix  which  satisfies  the  property  of  Eq.  (4.21),  i.e.  a  leading  matrix  whose  matrix  factors 
commute.  For  the  rigid-body  modes,  we  will  use  the  transition  matrix  already  developed  in  RA1. 
We  begin  by  denoting  the  complex  ( j  =  V-l)  eigenvalues  of  the  i*  mode  as 

=  -<;i©i +jcodi 

^7  =  _9i(0i  “j^di 

and  the  homogeneous  solution  to  Eq.  (4.14)  written  in  complex  form  is 

qj(t)  =  +  Cfe^1.  (4.29) 

where  C*  and  C~  are  constants  of  integration  for  the  ith  mode.  The  initial  condition  vector 
consistent  with  the  modal  state-vector  of  Eq.  (51)  is, 

*(.  =  0)  =  ^} 


and  hence  the  solution,  in  a  form  consistent  with  the  modal  state  vector  of  Eq.  (4.18),  becomes, 


qi(t)  = 


2co 


di 


eM(jq? +q?Ki  +  Ki^i))- Jextt(q? +q?(J®di  +?i®i)) 
*7eXit(jq?  +q?Ki  +j?i(0i))-  jAte^fq?  +q?(j(Odi  +qi®i)) 


(4.30) 


It  can  be  shown  that  the  modal  impulse  response  Hj  is  available  from  the  homogeneous 
solution  given  in  Eq.  (4.30)  using  the  following  initial  conditions: 


qi(t)  =  H;(t)  if 
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qi(t)  =  Hi(t)  if  q 


0_  J-tOi 
i_  1 


Due  to  its  importance  in  what  follows,  we  define  the  following  modal  state  initial  condition: 


We  now  define  the  following  matrices: 


[i-jSfe  -j— 

^di  ^di 

21+jSffi  j_L 

^di  ^di  _ 


Ai=diag(xt,>.T)  Pvi=AiPdi  (4.31) 


Given  these  definitions,  it  can  be  shown  that  the  modal  displacement  and  velocity  are. 


qi(t)  =  vTeAitPdiq? 
qi(t)  =  vTeAitPviq? 


(4.32) 


(4.33) 


where  vT  =  [1  1].  In  addition  to  the  relation  of  Eqs.  (4.31),  the  matrices  Pd  and  Pv  have  the 
additional  relation. 


Pvi  —  I*diAi 


(4.34) 


revealing  that  Aj  =  Aj  . 


Using  these  definitions,  we  can  write  the  total  modal  transient  response  as: 


qi(t)  =  vTeA^Pdiq?  +  j  vTeAi(t-'c)pdiqh^Tfi(T)dt 


(4.35) 


qj(t)  =  vTeAitPviq?  +  f  vTeAi(t 


(4.36) 
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It  is  of  interest  to  note  that  qh  =  B,  where  B  was  defined  in  RA1  above.  However,  it  should  be 
emphasized  that  q^  is  an  initial  condition  vector  which  produces  the  impulse  response  function 
from  the  homogeneous  solution  to  the  damped  modal  oscillator. 

We  now  define  a  complex  modal  state  vector  as  follows: 


q(0= 


qf(t) 


Using  this  state  vector,  Eqs.  (4.32)  and  (4.33)  can  be  written 


qi=vTqi 

qi=vT5i 


(4.37) 


(4.38) 

(4.39) 


where  we  have  dropped  the  time-dependence  notation  for  simplicity. 

We  now  construct  the  recursion  using  the  above  relations.  The  total  modal  response,  Eqs. 
(4.35)  and  (4.36)  can  be  written  in  terms  of  the  complex  modal  state  vector, 

<ii(t)  =  eA-‘Pdiq?  + 1  (4.40) 

0 

5i(t)  =  eAitPviq?+}eAi(t_T)Pviqh«hTf(t)dT  (4.41) 

0 

Equations  (4.40)  and  (4.41)  can  be  written  at  a  time  one  sample  At  later, 

t+At 

q1(t  +  at)  =  eA‘<t+i,)Pdiq?+  }  eA><,+4t-T)pdiqh*Tf(T)dt  (4.42) 

0 

t+At  , 

5i(t  +  At)  =  eA‘(t+At)Pviq?  +  j  eAi(t+At“I)pviqh<|)STf(x)dt  (4.43) 

0 

Using  the  fact  that  eAjt  satisfies  the  property  of  Eq.  (4.21),  Eqs.  (4.42)  and  (4.43)  become, 
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qi(t  +  At)  =  eAiAt 


:A'tPdiq?+jeA-('-'>Pd^hl|,Tf(T)dT 

0  J 


t+At 


...+ 


f  eA‘(l+4,“t)Pdiqh't|TfWdT 


(4.44) 


'  t 


4i  (t  +  At)  =  eA‘At  e AitPviqf  +  j  e Ai  ^  T)pviqh<fcTf (x)dx  + ... 

0  J 

t+At 

•••+  I  eAi(t+At“T)pviqh4>iTf(T)dx  (4.45) 

t 

If  we  again  assume  that  the  sampling  period  is  sufficiently  small  such  that  f  does  not  vary 
appreciably  over  At,  and  with  the  change  of  variable  c  =  t  +  At-x,  Eqs.  (4.44)  and  (4.45) 
become 

f  f  \ 


qi(t  +  At)  =  e 


A:  At 


eA|,Pdiq?+JeA‘(,-T)Pdiqh4|TfW*: 


0 


+  ... 


At 


) 


■ ...+  J  e^WdaPa^Ht) 


(4.46) 


qi(t  +  At)  =  e 


A:  At 


e  AitPviq?  +  j  e  Ai  (t~T)Pviqh<fcTf(x)dx 


At 


•••  +  J  eA‘HtoPvlqh*Tf(t) 
0 


(4.47). 


Recognizing  the  states  at  time  t  in  Eqs.  (4.46)  and  (4.47)  allows  the  following  simplification, 


At 


qi(t  +  At)  =  eA<A,qi(t)  +  J  eAi<<’Wpdiqh*lTf(t) 


0 

At 


5i(t  +  At)  =  eAiatqi(t)+  }  eAi('’>doPviqh*Tf(t) 

0 

and  the  integral  in  Eqs.  (4.48)  and  (4.49)  can  be  evaluated  as, 


(4.48) 

(4.49) 


At 


Sj(At)=  J  eA*(a)d<7  =  diag| 
0 


f  ^At_1  eXTAt_1' 


4  ’  ^7 
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which  yields  the  final  recursions  for  the  complex  modal  displacement  and  velocities 


qi(t  +  At)  =  Yiqi(t)  +  rdifi(t)  (4.50) 

(t  +  At)  =  (t)  +  T  vifj  (t)  (4.51) 

where  f{  =  <fcTf  and  in  which  the  following  quantities  are  defined: 

^i=eAiAt  (4.52) 

^\li  =  ^iPdiQh  (4.53) 

rvi  =  SjPyiQh  (4.54) 


The  explicit  form  of  RA2,  based  on  Eqs.  (4.25),  (4.50),  and  (4.51)  is  given  here.  For  clarity,  we 
again  indicate  the  various  quantities  at  time  kT  using  a  superscript  notation,  e.g.  q(kT)  =  qk.  The 
number  of  modes  retained  is  indicated  by  “p”  where  typically,  p«N. 


Recursive  Algorithm  #2  -  Explicit  Form 


Evaluate  Modal  Transition  Matrices  and  Input  Matrices  rd,  Tv 

Initialize  k  :=  0 
for  i  =  1 :  p 

?<k>:=*Tf(kT,x(k>',x<k>'l 
if i <  r  then 

q!k+1,  =  <Prbq!k)  +  rrbff° 


else 


end 


qi(k+1)  =  Tiqfk)  +  rdifp) 
^[k+1)  =  ^k)+rvif/k) 


end 

l 

X 


(k+1)  =  |  ^Vq(k+1) 

i=l  1 


>(k+l) 


T  T*(k+1) 


i=l 

k:=k+  1 
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4.4  Stability  Analysis  -  Recursive  Algorithm  1  -  Explicit  Form 


The  stability  analysis  for  the  explicit  form  of  RA1  is  given.  The  result  of  this  analysis 
reveals  a  fundamental  limitation  of  a  modal  transition  matrix  recursive  solution  to  Eq.  (2.19),  in 
that  the  spectral  radius  of  the  integration  matrix  operator  has  a  limiting  value  of  1.0,  approached 
from  above,  as  the  sample  length  is  made  smaller. 

For  the  stability  analysis,  the  structural  modifications  (e.g.  isolators)  are  assumed  to  be 
linear  elastic  and  hence  represented  by 

f(k)=-K*x*(k).  (4.55) 

The  linear  matrix  operator  representing  RA1  is  therefore, 

q(k+1>  =  [y-  r^TK*OTb]q(k)  =  Tiq(k)  (4.56) 

where  the  matrices  ¥  and  T,  defined  for  RA1,  are  of  size  2p  by  2p,  and  2p  by  p  respectively,  and 

contain  contributions  from  all  ‘p”  retained  modes,  and  Tb  (2p  by  2p)  is  a  Boolean  matrix  which 
rearranges  the  elements  of  q(k)  (2p  by  1)  to  be  consistent  with  form  of  Eq.  (4.55).  The  stability 
of  the  algorithm  is  measured  by  the  spectral  radius  of  T„  i.e.  stability  requires  that 

P(T,)  <  1. 

We  will  postpone  the  example  calculations  until  the  stability  analysis  of  RA2  is 
presented.  We  should  also  note  that  this  condition  on  the  spectral  radius  is  theoretical;  the  degree 
to  which  p  must  be  less  than  unity  is  dependent  on  numerical  considerations. 


4.5  Stability  Analysis  -  Recursive  Algorithm  2  -  Explicit  Form 

We  define  a  combined  complex  modal  state  vector,  of  dimensions  4x1,  as 

«“■(* 


(4.57) 
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and  along  with  Eq.  (4.55),  allows  the  elastic  mode  recursion  of  RA2  to  be  written  as  the 
following  linear  matrix  operator, 

Q(k+1)  _  _  roTK*<l>VT]Q(k)  =  T2Q(k)  (4.58) 

where  the  matrices  Y  and  T,  defined  for  RA2,  are  4p  by  4p,  and  4p  by  p,  respectively,  and  VT  is 

p  by  4p.  Each  of  these  matrices  is  block  quasi-diagonal  where  each  block  is  associated  with  a 
retained  mode,  and  the  definition  of  which  can  be  found  in  the  above  section  outlining  the 
development  of  RA2.  The  stability  of  the  algorithm  is  measured  by  the  spectral  radius  of  T2,  i.e. 
stability  requires  that  p(T2)  <  1. 

4.6  Example  Calculations  -  Stability  of  Explicit  Recursive  Algorithms 

We  will  make  use  of  a  simple  4  DOF  spring-mass  model  to  calculate  algorithm  stability. 
The  system  is  shown  in  Figure  9.  In  this  example,  the  synthesis  will  install  the  stiffness 
modification  k*,  and  calculate  the  transient  response  of  mass  #1,  mi-  The  example  system  has  the 
following  parameter  values, 

M  =  diag(1.7513  2.1016  0.3503  1.7513)  kg, 

k*  =  700  N/m  k,  =  k2  =  k3  =  175. 13  N/m, 

and  has  natural  frequencies  in  Hertz,  f3  =  0.6261,  f2=  1.5421,  f3  =  2.5750  f4=  5.2678  (Hz), 
k* 


^wv- 


Figure  9.  A  Lumped  4-DOF  System  with  a  stiffness  modification 

As  the  purpose  of  this  example  is  to  determine  stability  of  the  algorithms  RA1  and  RA2,  we  limit 
the  example  to  a  linear  elastic  modification. 
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Stability  of  RA1:  We  first  calculate  the  spectral  radius  of  the  explicit  operator  for  RA1  (Eq. 

(4.56)),  as  a  function  of  sample  length.  At,  for  a  fixed  value  of  k*  =  700  N/m,  with  C,  =  0.0, 

C  =  0.05,  C  =  0.2,  and  £  =  0.3.  This  is  shown  in  Figure  10.  The  range  of  sample  lengths  is  from 

0.0001  s  to  0.01  s,  which  is  consistent  with  the  range  of  natural  periods  for  the  model  which  is 
0.19  s  to  1.6  s,  with  respect  to  the  requirements  on  a  sufficiently  small  time  step. 


Figure  10.  Spectral  radius  of  Tj  versus  sample  length  At  for  various  £ 

From  this  plot,  it  is  clear  that  the  simple  modal-based  recursion  is  ineffective  for  the  solution  of 
Eq.  (2.14).  The  algorithm  achieves  marginal  stability,  as  the  step  size  is  made  small,  for  large 
values  of  modal  damping.  This  result  is  attributed  to  the  use  of  the  modal  transition  matrix,  eAAt, 
which  clearly  has  a  limiting  spectral  radius  of  1,  for  At  0.  By  increasing  the  modal  damping, 

the  region  of  neutral  stability  is  extended  to  longer  sample  lengths,  but  this  does  not  offer  any 
practical  value. 

We  can  also  examine  the  dependence  of  the  stability  on  the  stiffness  value  k*,  for  fixed 
sample  length  of  At  =  0.01s,  and  for  various  values  of  £,  shown  in  Figure  11.  Again, 

unrealistically  high  values  of  modal  damping  are  required  to  stabilize  the  algorithm. 


61 


1.015 
1.01 
1.005 

1 

0.995 
0.99 
0.985 

0  200  400  600  800 

Stiffness  k*  (N/m)  with  At  =  0.01  sec 

Figure  11.  Spectral  radius  of  T,  versus  stiffness  k*  for  various  C, 

Stability  of  RA2:  We  calculate  the  spectral  radius  of  the  operator  T,  (Eq.  (4.58))  for 

explicit-RA2,  as  a  function  of  sample  length,  At,  for  k*  =  700  N/m,  and  for  various  values  of  C,. 

The  results  of  this  calculation  are  shown  in  Figure  12.  These  results  are  in  fact  identical  to  the 
corresponding  results  for  RA1,  indicating  that  the  use  of  the  modal  transition  matrix  based  on  the 
second-order  equations  does  not  have  any  effect  on  the  stability,  as  might  be  expected,  as  the 
difference  between  the  algorithms  is  essentially  a  change  of  coordinates,  and  both  algorithms  are 
explicit.  It  is  seen  that  RA2  cannot  be  stabilized  by  increasing  the  modal  damping. 

The  dependence  of  stability  of  RA2  with  varying  k*  is  shown  in  Figure  13,  and  again  the 
results  of  Figure  13  are  identical  to  those  calculated  using  RA1,  while  varying  k*. 
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Figure  12.  Spectral  radius  of  T2  versus  sample  length  At  for  various  £ 


Figure  13.  Spectral  radius  of  T2  versus  stiffness  k*  for  various  C, 

■> 

The  algorithms  RA1  and  RA2  are  attractive  in  that  they  satisfy  some  of  the  requirements 
set  forth  above  regarding  the  development  of  a  solution  method  for  Eq.  (2.14),  specifically  in  that 
these  algorithms  achieve  a  model  reduction  using  modal  superposition.  The  results  shown 
immediately  above  indicate  that  these  algorithms  are  not  viable  methods,  in  their  explicit  forms. 
We  will  now  develop  implicit,  predictor-corrector  versions  of  these  algorithms,  to  be  referred  to 
as  RA1PC,  and  RA2PC,  to  see  if  the  instability  of  the  explicit  forms  can  be  overcome. 
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4.7  Recursive  Algorithm  1:  Predictor-Corrector  Form 


We  begin  with  a  recursive  form  of  the  solution  to  Eq.  (4.14) 

q|k+1)  =  eAiAtq|k>  +  J  eA^t+At_T^Bfidx, 

t 

and  letting  a  =  t  +  At  -  x,  then  da  =  -dx  and  Eq.  (4.59)  becomes, 

q[k+1)  =  eAi Atq[k)  +  j  eAi<TBfj  (t  +  At  -  x)da . 

0 

We  employ  the  following  interpolation  for  the  modal  force, 

fi  =  afj  (t)  +  (1  -  oc)fj  (t  +  At) 
where  0  <  a  <  1.  Using  a  =  a/At,  we  have, 

q-k+1)  =eAiAtq[k)  +  J  eA>cB|-^[^af^k)  +(At-a)f^k+1)]|da, 


or, 


n(k+l)  =  eAiAtn(k>  +  —I 


At  .  )  ,  ...  _..  (At  , 

J  eAi(Tada  B(f^ -f^k+1U  +  J  eAi<Tda 
lo  J  v  '  VO 


+ 


At 


Bf<k+1>. 


The  integrals  in  Eq.  (4.63)  can  now  be  evaluated. 


At 

\ 

0 

At 

f 
0 

At 

J 
0 


JeAiado  =  Af1(eAiAt-l), 

0  \J 

Tejada  =  AtAjV1*  -  Aj-2(eA‘At  - 1). 

o  v  ' 

Letting  *P  =  eA;At  and  T  =  J  eAi°da  =  Af^eAiAt  -  ijfi,  then. 


(4.59) 


(4.60) 


(4.61) 


(4.62) 


(4.63) 


q-k+1)  =  ^qSk)  +r^k+1)  +  A"1YiB(f^k)  -f^^j-^A"1^^  -f<k+1))  (4.64) 


or  finally. 
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qjk+1)  =  ^q5k)  +  - — TjVf^  -  f<k+1))  +  r$k+1) 


(4.65) 


The  predictor-corrector  algorithm,  RA1PC  is  now  outlined. 


Recursive  Algorithm  #1  -  Predictor-Corrector 


•  k  :=  0 


Predictor: 


Corrector: 


for  i  =  1  :  p 

q-k+l)P  =YiqJk)+rifp) 


£(k+1)p  _  ^q(k+l)p 


x(k+l)p  =  <&q(k+1)P 
for  i  =  1 :  p 


f/k+1)'=/x(k+1)P,x(k+1)' 


for  i  =  1 :  p 


q[k+1)  =  Yiqjk)  +  TjYfp)  -f<k+1)P)  +rjf/k+1) 


x(k+1)  _  ^q(k+l) 
x(k+1)  _  ^q(k+l) 

for  i  =  1 :  p 


f<k+1)=f(x(k+1),x(k+1) 


end 

k  :=  k  +  1 
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4.8  Stability  Analysis  -  Recursive  Algorithm  1  Predictor-Corrector  Form 

For  the  linear  stability  analysis,  we  will  again  assume  that 

f(k)  =  -K*x*(k).  (4.66) 

Combining  the  equations  above  yields  the  following, 

q(k+1)=[(I-rZ)T+(rZ)2+0Z(^-rZ-I)]q(k)  (4.67) 

T 

where  Z  =  <I>K  *  <J>Tq ,  and  Tq  is  a  Boolean  matrix  which  extracts  displacement  entries  in  the  q 
vector.  Equation  (4.67)  can  be  abbreviated  as 

q[k+1>  =T1pcqJk),  (4.68) 

and  the  condition  for  stability  is  that  the  spectral  radius  of  Tjpc  does  not  exceed  unity,  i.e. 
P(T1pc)<l. 


4.9  Recursive  Algorithm  2:  Predictor-Corrector  Form 


We  begin  with  the  following  equations  for  the  complex  modal  state  responses, 


qp+1>  =  e'^qp)  +  J  eAi(t+At-T)pdiqhfj(x)dT 


(4.69) 


^0  =  eA.^k)+  |  eA,(.+At-t)Pv.-h?.(t)(h. 


(4.70) 


Making  use  of  Eq.  (4.61),  and  following  the  same  procedure  as  used  in  RA1PC,  we  arrive  at  the 
following  recursions  for  the  i*  complex  modal  state  responses. 


qp+1) = 'F,5ik)  +  edl(  f»>  -  fp+»p) + rdif[k- 


(4.71) 


-(k+1)  =  qr-p)  +  q( fp)  _  fp+l)p"j  +  rvif^k+1)^ 


(4.72) 


where  the  following  quantities  are  defined: 

%  =  eAiAt 


(4.73) 
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rdi=AT1(>ri-I)Pdiqh  (4.74) 

rv,  =AT1(4',-I)Pviq1,  (4.75) 

e<ii=AT1^iPdiqh-.lrdij  (4.76) 

Syi  -  A7*^'FiPviqh  -  ■^■rvi j  (4.77) 


The  predictor-corrector  algorithm,  RA2PC  is  now  outlined.  In  this  algorithm,  the  relevant 
matrices  have  been  expanded  to  include  all  retained  modes;  that  is,  we  eliminate  the  “for-end” 
loops  over  the  modes  in  favor  of  a  matrix  notation,  where  the  various  matrices  are  block- 
diagonal. 


Recursive  Algorithm  #2  -  Predictor-Corrector 


•  k  :=  0 
Predictor: 

•  Q(k+1)P  =  'J'qM  +  rdf  ^k+1^ 

.  <J(k+1)P  =  +  rdif(k+1) 

.  x(k+l)P=oyTQ(k+l)p 

x(k+1)P  =  OVvTQ(k+1>P 

•  for  i  =  1 :  p 

f<k+1>P=ffx<k+1>P,x<k+1>P) 

•  end 
Corrector: 

•  Q(k+1)  =  *PQ(k)  +  0d^f  ^  -  f 1 (k+1)P  j  +  Fdf(k+1)P 

.  $(k+1)  =  ^$(k) + ev(^f  (k>  -  f  (k+1)p  j  +  rvf {k+1)p 

.  x(k+1)p=4>vdTQ(k+1)P 

x(k+1>P  =  <&vjQ(k+1)P 

•  for  i  =  1 :  p 

f<k+»  =  f(x<k+1>,x<k+1>) 

•  end 

•  k  :=  k  +  1 
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4.10  Stability  Analysis  -  Recursive  Algorithm  2  Predictor-Corrector  Form 

For  the  linear  stability  analysis,  we  will  assume  that 

f(k)=-K*x*(k)  (4.78) 

Combining  the  equations  above  yields  the  following, 

Q(k+1) = [(i  -  rz)*F + (rz)2 + ©z(t  -  rz  -  i)]Q(k)  (4.79) 

T  T 

where  Z  =  <&  K  *  4>V  ,  and  V  is  a  Boolean  matrix  which  effects  the  addition  of  the  elements  of 
the  complex  modal  responses.  Equation  (4.79)  can  be  written  as 

Q(k+1)  _  TPCQ(k)  (4.80) 

The  theoretical  condition  for  stability  is  that  the  spectral  radius  of  T|*c  does  not  exceed  unity,  i.e 

P(T2pc)<l. 

4.11  Example  Calculations  -  Stability  of  Predictor-Corrector  Algorithms 

The  calculations  will  make  use  of  the  system  of  Figure  9,  with  identical  parameters  for  this 
example. 

Stability  of  RA1PC:  We  first  calculate  the  spectral  radius  of  the  RA1PC  operator,  (Eq.  (4.68)), 
as  a  function  of  sample  length.  At,  for  a  fixed  value  of  k*  =  700  N/m,  with  £  =  0.0,  C,  =  0.05,  C,  = 

0.2,  and  C,  =  0.3.  This  is  shown  in  Figure  14.  The  range  of  sample  lengths  is  from  0.0001  s  to 
0.01  sec.  The  range  of  natural  periods  for  the  model  is  0.19  s  to  1.6  s. 
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Figure  14.  Spectral  radius  of  versus  sample  length  At  for  various  £ 

From  Figure  14,  it  is  seen  that,  as  compared  with  the  explicit  algorithm  RA1  (Figure  10),  the 
reformulation  of  the  basic  algorithm  as  an  implicit  predictor-corrector,  using  a  single  corrector 
step,  has  provided  a  substantial  improvement  in  the  fundamental  stability  of  the  algorithm.  The 
spectral  radius  decreases  with  increasing  step  size,  for  non-zero  system  modal  damping.  The 
dependence  of  stability  of  RA1PC  with  varying  k*  is  shown  in  Figure  15,  for  At  =  0.01s. 


Figure  15.  Spectral  radius  of  T^0  versus  sample  length  K*  for  various  C, 
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Stability  of  RA2PC:  We  calculate  the  spectral  radius  of  the  operator  for  RA2PC  as  a  function  of 
sample  length,  At,  for  a  fixed  value  of  k*  =  700  N/m,  with  £  =  0.0,  £  =  0.05,  C,  =  0.2,  and  £  =  0.3. 

This  is  shown  in  Figure  16.  The  range  of  sample  lengths  is  from  0.0001  s  to  0.01  sec.  The  range 
of  natural  periods  for  the  model  is  0.19  s  to  1.6  s.  The  reformulation  of  RA2  as  an  implicit 
algorithm  has  again  provided  a  substantial  improvement  in  the  stability.  As  seen  from  Figure  16, 
the  calculation  produces  results  identical  to  those  from  RA1PC,  again  revealing  the  fundamental 
similarity  of  the  algorithms,  the  difference  being  the  use  of  the  complex  modal  state  vector  in 
RA2  versus  the  traditional  modal  state  vector  in  RA1.  The  dependence  of  stability  of  RA2PC 
with  varying  k*  is  shown  in  Figure  17,  for  At  =  0.01s,  which  is  again,  identical  to  the  analogous 

calculation  for  RA1PC. 


Figure  16.  Spectral  radius  of  T|c  versus  sample  length  At  for  various  C, 
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Figure  17.  Spectral  radius  of  X£c  versus  sample  length  k*  for  various 


Spectral  Radius  (sec'2)  vs  k* 


4. 12  Transient  Synthesis  Example  -  Linear  Elastic  Modification 


The  algorithms  RA1PC  and  RA2PC  will  be  used  to  synthesize  the  transient  response  of 
the  system  shown  in  Figure  18.  The  synthesis  will  be  used  to  simultaneously  replace  spring  kj 
with  k*  and  calculate  the  response  of  mass  m,  to  the  ground  displacement  input  y(t).  The  ground 
motion  y(t)  is  taken  as  a  unit  step  input. 


k* 


-JVW- 


HWH mi  H\MH m2 1 — VW — m3  ~ ^WV-  m4 

kl  k2 


k3 


k4 


->y(t)  =  l(t>0) 


Figure  18.  System  for  synthesis  to  transient  base  displacement  response  with  spring  replacement 


The  purpose  of  this  example  is  two-fold.  The  first  purpose  is  to  demonstrate  that  a  single 
corrector  step  taken  in  algorithms  RA1PC  and  RA2PC,  while  providing  a  spectral  radius  less 
than  unity,  does  not  provide  a  sufficiently  small  spectral  radius  to  stabilize  the  algorithms  for  a 
range  of  stiffness  modification  magnitudes.  The  recursive  algorithm  solutions  will  be  compared 
with  a  direct  integration  (variable  time  step  Runge-Kutte)  of  the  system  equations  for  the 


71 


modified  system,  i.e.  lq  =  k*.  In  these  calculations,  the  synthesis  is  first  performed  using  a  value 
of  k*  =  k,  =  k3  =  k4  =  175.13  N/m,  with  £  =  0.0.  The  results  from  RA1PC  are  shown  in  Figure  19, 

and  from  RA2PC  in  Figure  20,  where  both  algorithms  provide  accurate  (and  identical)  results. 
The  next  calculations  use  k*  =  1.5k,  =  262.70  N/m,  and  the  instability  is  evident  in  the  results 
from  both  algorithms,  shown  in  Figures  21  and  22.  These  calculations  with  k*  =  1.5k,  are  then 
repeated  using  £  =  0.05,  and  the  system  modal  damping  stabilizes  the  algorithm,  as  shown  in 

Figures  23  and  24. 

The  second  purpose  of  this  example  is  to  demonstrate  the  substantial  decrease  in 
computer  time  required  by  the  algorithms  based  on  the  complex  modal  state  vector 
(RA2/RA2PC).  Here,  algorithms  RA1  and  RA2  will  be  compared  to  a  standard  modal  approach 
to  the  solution  of  this  system  with  k*  already  installed,  as  well  as  to  a  direct  integration  of  the 
equations  of  motion.  In  other  words,  we  will  be  comparing  the  computer  time  required  for  a 
structural  modification  calculation  with  that  required  for  a  standard  analysis,  where  the  structural 
modification  has  already  been  installed  during  the  model  assembly  phase,  i.e.  kt  =  k*.  As  the 
intended  use  of  these  algorithms  is  for  the  solution  of  nonlinear  structural  synthesis  problems,  the 
comparison  includes  the  direct  integration  solution  time.  The  modal  solution  is  included  as  a 
reference  time  for  a  standard  linear  solution.  The  analysis  corresponds  to  that  analysis  whose 
results  are  shown  in  Figure  24,  but  using  the  explicit  form,  RA1  and  RA2.  From  the  timing 
results  from  Table  3,  it  is  seen  that  the  time  required  for  the  synthesis  algorithm  RA2  is  of  the 
same  order  of  magnitude  as  the  standard  linear  modal  transient  analysis.  The  algorithm  RA1 
requires  a  time  comparable  to  the  direct  integration.  Of  course,  as  the  model  size  increases,  the 
direct  integration  time  increases,  while  the  times  for  RA1  and  RA2  are  independent  of  model 
size,  once  a  modal  database  has  been  calculated. 


Table  3.  Comparison  of  Times  Required  for  Algorithms 


Algorithm: 

Direct  Integ. 

Modal 

RA1 

RA2 

Time  (sec) 

1.130 

1.019 

0.175 
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Stability  Comparison:  Transient  Base  Displacement  Response  At 
Direct  Integration -  Synthesis  — 


lUft 


Time  (sec) 


Time  (sec) 


Figure  19.  RA1PC:  k*  =  175. 13  N/m,  £  =  0.0  Figure  20.  RA2PC:  k*  =  175. 13  N/m,  £  =  0.0 


2  3 

Time  (sec) 


Time  (sec) 


Figure  21 .  RA1PC:  k*  =  262.70  N/m,  £  =  0.0  Figure  22.  RA2PC:  k*  =  262.70  N/m,  C  =  0.0 


aA 


A 

W* 


2  3 

Time  (sec) 


2  3 

Time  (sec) 


Figure  23.  RA1PC:  k*  =  262.70  N/m,  £  =  0.05  Figure  24.  RA2PC:  k*  =  262.70  N/m,  C  =  0.05 


4. 13  Summary  of  Results  for  Transition  Matrix-Based  Recursive  Algorithms: 
Explicit  and  Predictor-Corrector 

The  goal  of  this  work  is  the  development  of  an  efficient  and  robust  algorithm  for  the 
solution  of  the  governing  integral  equation  for  locally  nonlinear  transient  structural  synthesis. 
The  new  algorithm  must  preserve  the  existing  features  of  the  existing  formulation.  The  features 
include, 

•  a  solution  time  which  is  independent  of  model  size;  an  exact  model  reduction  which  is 
essentially  unrestricted  in  the  dimension  of  the  reduction  is  implicit  in  the  formulation, 

•  the  ability  to  perform  local  structural  modifications  of  general  nonlinear  characteristics 

•  accurate  and  fast  solutions  facilitating  optimal  design 

Two  general  algorithms  have  developed  and  analyzed  for  stability.  Both  algorithms  are 
based  on  modal  transition  matrices,  facilitating  the  recursive  forms  of  these  algorithms.  The  first 
algorithm,  referred  to  as  Recursive  Algorithm  #1  (RA1),  is  based  on  the  standard  modal  state 
vector  and  the  associated  transition  matrix,  e^1.  The  second  algorithm  is  based  on  a  newly 

defined  complex  modal  state  vector,  and  also  defines  a  new  transition  matrix  formulation  based 
exclusively  on  the  second-order  modal  differential  equations.  This  algorithm  is  referred  to  as 
Recursive  Algorithm  #2  (RA2).  It  was  shown  that  while  the  formulations  of  these  algorithms  are 
different,  RA1  being  a  real-arithmetic  algorithm  and  RA2  being  complex,  the  algorithms  produce 
identical  numerical  results.  However,  RA2  allows  a  greater  level  of  “diagonalization”  relative  to 
RA1,  and  hence  yields  computing  times  of  an  order  of  magnitude  less  than  that  required  by  RA1 
(see  Table  1). 

The  inherent  instability  of  algorithms  RA1  and  RA2  was  addressed  by  their  reformulation 
as  predictor-corrector  algorithms,  using  a  linear  interpolation  for  the  modal  force,  and  a  single 
corrector  step.  This  reformulation  was  successful  in  that  it  demonstrated  that  the  limiting  value 
of  unity  for  the  spectral  radii  of  RA1  and  RA2  could  be  removed,  and  that  spectral  radii  less  than 
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unity  are  achievable.  However,  the  study  performed  indicated  that  the  degree  of  stability 
provided  by  this  basic  reformulation  is  not  sufficient  to  provide  a  robust  algorithm.  Therefore, 
rather  than  continue  with  the  development  of  these  recursive  predictor-corrector  algorithms,  this 
work  will  continue  with  the  development  of  a  recursive  algorithm  which  is  not  based  on  a 
transition  matrix,  but  rather  on  a  block-by-block  convolution  algorithm.  It  will  be  shown  in  the 
next  section  that  this  block-by-block  convolution  algorithm  preserves  all  the  features  of  the 
existing  synthesis  formulation,  and  provides  a  dramatic  decrease  in  the  computing  time  required. 
This  algorithm  is  considered  a  non-modal  formulation,  as  it  is  based  on  the  use  of  physical 
impulse  response  to  represent  the  (linear)  substructures,  which  are  subjected  to  structural 
modification,  substructure  coupling,  and  ground-motion  input  through  isolation.  This  algorithm 
will  be  shown  to  be  exponentially  convergent,  regardless  of  the  linear  substructure  properties, 
and  for  a  general  and  broad  class  of  nonlinear  modifications. 
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